Numerical Analysis

• 2.I.5E

  Numerical Analysis

  Part IB, 2001

  Find an LU factorization of the matrix

  $$A=\left(\begin{array}{rrrr} 2 & -1 & 3 & 2 \\ -4 & 3 & -4 & -2 \\ 4 & -2 & 3 & 6 \\ -6 & 5 & -8 & 1 \end{array}\right)$$

  and use it to solve the linear system $A \mathbf{x}=\mathbf{b}$, where

  $$\mathbf{b}=\left(\begin{array}{r} -2 \\ 2 \\ 4 \\ 11 \end{array}\right)$$

  comment

• 2.II.14E

  Numerical Analysis

  Part IB, 2001

  (a) Let $B$ be an $n \times n$ positive-definite, symmetric matrix. Define the Cholesky factorization of $B$ and prove that it is unique.

  (b) Let $A$ be an $m \times n$ matrix, $m \geqslant n$, such that $\operatorname{rank} A=n$. Prove the uniqueness of the "skinny QR factorization"

  $$A=Q R,$$

  where the matrix $Q$ is $m \times n$ with orthonormal columns, while $R$ is an $n \times n$ upper-triangular matrix with positive diagonal elements.

  [Hint: Show that you may choose $R$ as a matrix that features in the Cholesky factorization of $B=A^{T} A$.]

  comment

• 3.I.6E

  Numerical Analysis

  Part IB, 2001

  Given $f \in C^{n+1}[a, b]$, let the $n$ th-degree polynomial $p$ interpolate the values $f\left(x_{i}\right)$, $i=0,1, \ldots, n$, where $x_{0}, x_{1}, \ldots, x_{n} \in[a, b]$ are distinct. Given $x \in[a, b]$, find the error $f(x)-p(x)$ in terms of a derivative of $f$.

  comment

• 3.II.16E

  Numerical Analysis

  Part IB, 2001

  Let the monic polynomials $p_{n}, n \geqslant 0$, be orthogonal with respect to the weight function $w(x)>0, a<x<b$, where the degree of each $p_{n}$ is exactly $n$.

  (a) Prove that each $p_{n}, n \geqslant 1$, has $n$ distinct zeros in the interval $(a, b)$.

  (b) Suppose that the $p_{n}$ satisfy the three-term recurrence relation

  $$p_{n}(x)=\left(x-a_{n}\right) p_{n-1}(x)-b_{n}^{2} p_{n-2}(x), \quad n \geqslant 2$$

  where $p_{0}(x) \equiv 1, p_{1}(x)=x-a_{1}$. Set

  $$A_{n}=\left(\begin{array}{ccccc} a_{1} & b_{2} & 0 & \cdots & 0 \\ b_{2} & a_{2} & b_{3} & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & b_{n-1} & a_{n-1} & b_{n} \\ 0 & \cdots & 0 & b_{n} & a_{n} \end{array}\right), \quad n \geqslant 2 .$$

  Prove that $p_{n}(x)=\operatorname{det}\left(x I-A_{n}\right), n \geqslant 2$, and deduce that all the eigenvalues of $A_{n}$ reside in $(a, b)$.

  comment

• 2.I.5B

  Numerical Analysis

  Part IB, 2002

  Applying the Gram-Schmidt orthogonalization, compute a "skinny"

  QR-factorization of the matrix

  $$A=\left[\begin{array}{lll} 1 & 1 & 2 \\ 1 & 3 & 6 \\ 1 & 1 & 0 \\ 1 & 3 & 4 \end{array}\right],$$

  i.e. find a $4 \times 3$ matrix $Q$ with orthonormal columns and an upper triangular $3 \times 3$ matrix $R$ such that $A=Q R$.

  comment

• 2.II.14B

  Numerical Analysis

  Part IB, 2002

  Let $f \in C[a, b]$ and let $n+1$ distinct points $x_{0}, \ldots, x_{n} \in[a, b]$ be given.

  (a) Define the divided difference $f\left[x_{0}, \ldots, x_{n}\right]$ of order $n$ in terms of interpolating polynomials. Prove that it is a symmetric function of the variables $x_{i}, i=0, \ldots, n$.

  (b) Prove the recurrence relation

  $$f\left[x_{0}, \ldots, x_{n}\right]=\frac{f\left[x_{1}, \ldots, x_{n}\right]-f\left[x_{0}, \ldots, x_{n-1}\right]}{x_{n}-x_{0}}$$

  (c) Hence or otherwise deduce that, for any $i \neq j$, we have

  $f\left[x_{0}, \ldots, x_{n}\right]=\frac{f\left[x_{0}, \ldots, x_{i-1}, x_{i+1}, \ldots, x_{n}\right]-f\left[x_{0}, \ldots, x_{j-1}, x_{j+1}, \ldots, x_{n}\right]}{x_{j}-x_{i}} .$

  (d) From the formulas above, show that, for any $i=1, \ldots, n-1$,

  $$f\left[x_{0}, \ldots, x_{i-1}, x_{i+1}, \ldots, x_{n}\right]=\gamma f\left[x_{0}, \ldots, x_{n-1}\right]+(1-\gamma) f\left[x_{1}, \ldots, x_{n}\right],$$

  where $\gamma=\frac{x_{i}-x_{0}}{x_{n}-x_{0}}$.

  comment

• 3.I.6B

  Numerical Analysis

  Part IB, 2002

  For numerical integration, a quadrature formula

  $$\int_{a}^{b} f(x) d x \approx \sum_{i=0}^{n} a_{i} f\left(x_{i}\right)$$

  is applied which is exact on $\mathcal{P}_{n}$, i.e., for all polynomials of degree $n$.

  Prove that such a formula is exact for all $f \in \mathcal{P}_{2 n+1}$ if and only if $x_{i}, i=0, \ldots, n$, are the zeros of an orthogonal polynomial $p_{n+1} \in \mathcal{P}_{n+1}$ which satisfies $\int_{a}^{b} p_{n+1}(x) r(x) d x=0$ for all $r \in \mathcal{P}_{n}$. [You may assume that $p_{n+1}$ has $(n+1)$ distinct zeros.]

  comment

• 3.II.16B

  Numerical Analysis

  Part IB, 2002

  (a) Consider a system of linear equations $A x=b$ with a non-singular square $n \times n$ matrix $A$. To determine its solution $x=x^{*}$ we apply the iterative method

  $$x^{k+1}=H x^{k}+v .$$

  Here $v \in \mathbb{R}^{n}$, while the matrix $H \in \mathbb{R}^{n \times n}$ is such that $x^{*}=H x^{*}+v$ implies $A x^{*}=b$. The initial vector $x^{0} \in \mathbb{R}^{n}$ is arbitrary. Prove that, if the matrix $H$ possesses $n$ linearly independent eigenvectors $w_{1}, \ldots, w_{n}$ whose corresponding eigenvalues $\lambda_{1}, \ldots, \lambda_{n}$ satisfy $\max _{i}\left|\lambda_{i}\right|<1$, then the method converges for any choice of $x^{0}$, i.e. $x^{k} \rightarrow x^{*}$ as $k \rightarrow \infty$.

  (b) Describe the Jacobi iteration method for solving $A x=b$. Show directly from the definition of the method that, if the matrix $A$ is strictly diagonally dominant by rows, i.e.

  $$\left|a_{i i}\right|^{-1} \sum_{j=1, j \neq i}^{n}\left|a_{i j}\right| \leq \gamma<1, \quad i=1, \ldots, n,$$

  then the method converges.

  comment

• 2.I.5B

  Numerical Analysis

  Part IB, 2003

  Let

  $$A=\left(\begin{array}{cccc} 1 & a & a^{2} & a^{3} \\ a^{3} & 1 & a & a^{2} \\ a^{2} & a^{3} & 1 & a \\ a & a^{2} & a^{3} & 1 \end{array}\right), \quad b=\left(\begin{array}{c} \gamma \\ 0 \\ 0 \\ \gamma a \end{array}\right), \quad \gamma=1-a^{4} \neq 0$$

  Find the LU factorization of the matrix $A$ and use it to solve the system $A x=b$.

  comment

• 2.II.14B

  Numerical Analysis

  Part IB, 2003

  Let

  $$f^{\prime \prime}(0) \approx a_{0} f(-1)+a_{1} f(0)+a_{2} f(1)=\mu(f)$$

  be an approximation of the second derivative which is exact for $f \in \mathcal{P}_{2}$, the set of polynomials of degree $\leq 2$, and let

  $$e(f)=f^{\prime \prime}(0)-\mu(f)$$

  be its error.

  (a) Determine the coefficients $a_{0}, a_{1}, a_{2}$.

  (b) Using the Peano kernel theorem prove that, for $f \in C^{3}[-1,1]$, the set of threetimes continuously differentiable functions, the error satisfies the inequality

  $$|e(f)| \leq \frac{1}{3} \max _{x \in[-1,1]}\left|f^{\prime \prime \prime}(x)\right| .$$

  comment

• 3.I.6B

  Numerical Analysis

  Part IB, 2003

  Given $(n+1)$ distinct points $x_{0}, x_{1}, \ldots, x_{n}$, let

  $$\ell_{i}(x)=\prod_{\substack{k=0 \\ k \neq i}}^{n} \frac{x-x_{k}}{x_{i}-x_{k}}$$

  be the fundamental Lagrange polynomials of degree $n$, let

  $$\omega(x)=\prod_{i=0}^{n}\left(x-x_{i}\right)$$

  and let $p$ be any polynomial of degree $\leq n$.

  (a) Prove that $\sum_{i=0}^{n} p\left(x_{i}\right) \ell_{i}(x) \equiv p(x)$.

  (b) Hence or otherwise derive the formula

  $$\frac{p(x)}{\omega(x)}=\sum_{i=0}^{n} \frac{A_{i}}{x-x_{i}}, \quad A_{i}=\frac{p\left(x_{i}\right)}{\omega^{\prime}\left(x_{i}\right)}$$

  which is the decomposition of $p(x) / \omega(x)$ into partial fractions.

  comment

• 3.II.16B

  Numerical Analysis

  Part IB, 2003

  The functions $H_{0}, H_{1}, \ldots$ are generated by the Rodrigues formula:

  $$H_{n}(x)=(-1)^{n} e^{x^{2}} \frac{d^{n}}{d x^{n}} e^{-x^{2}}$$

  (a) Show that $H_{n}$ is a polynomial of degree $n$, and that the $H_{n}$ are orthogonal with respect to the scalar product

  $$(f, g)=\int_{-\infty}^{\infty} f(x) g(x) e^{-x^{2}} d x$$

  (b) By induction or otherwise, prove that the $H_{n}$ satisfy the three-term recurrence relation

  $$H_{n+1}(x)=2 x H_{n}(x)-2 n H_{n-1}(x) .$$

  [Hint: you may need to prove the equality $H_{n}^{\prime}(x)=2 n H_{n-1}(x)$ as well.]

  comment

• 2.I.9A

  Numerical Analysis

  Part IB, 2004

  Determine the coefficients of Gaussian quadrature for the evaluation of the integral

  $$\int_{0}^{1} f(x) x d x$$

  that uses two function evaluations.

  comment

• 2.II.20A

  Numerical Analysis

  Part IB, 2004

  Given an $m \times n$ matrix $A$ and $\mathbf{b} \in \mathbb{R}^{m}$, prove that the vector $\mathbf{x} \in \mathbb{R}^{n}$ is the solution of the least-squares problem for $A \mathbf{x} \approx \mathbf{b}$ if and only if $A^{T}(A \mathbf{x}-\mathbf{b})=\mathbf{0}$. Let

  $$A=\left[\begin{array}{cc} 1 & 2 \\ -3 & 1 \\ 1 & 3 \\ 4 & 1 \end{array}\right], \quad \mathbf{b}=\left[\begin{array}{c} 3 \\ 0 \\ -1 \\ 2 \end{array}\right]$$

  Determine the solution of the least-squares problem for $A \mathbf{x} \approx \mathbf{b}$.

  comment

• 3.I.11A

  Numerical Analysis

  Part IB, 2004

  The linear system

  $$\left[\begin{array}{lll} \alpha & 2 & 1 \\ 1 & \alpha & 2 \\ 2 & 1 & \alpha \end{array}\right] \mathbf{x}=\mathbf{b}$$

  where real $\alpha \neq 0$ and $\mathbf{b} \in \mathbb{R}^{3}$ are given, is solved by the iterative procedure

  $$\mathbf{x}^{(k+1)}=-\frac{1}{\alpha}\left[\begin{array}{lll} 0 & 2 & 1 \\ 1 & 0 & 2 \\ 2 & 1 & 0 \end{array}\right] \mathbf{x}^{(k)}+\frac{1}{\alpha} \mathbf{b}, \quad k \geqslant 0$$

  Determine the conditions on $\alpha$ that guarantee convergence.

  comment

• 3.II.22A

  Numerical Analysis

  Part IB, 2004

  Given $f \in C^{3}[0,1]$, we approximate $f^{\prime}\left(\frac{1}{3}\right)$ by the linear combination

  $$\mathcal{T}[f]=-\frac{5}{3} f(0)+\frac{4}{3} f\left(\frac{1}{2}\right)+\frac{1}{3} f(1)$$

  By finding the Peano kernel, determine the least constant $c$ such that

  $$\left|\mathcal{T}[f]-f^{\prime}\left(\frac{1}{3}\right)\right| \leq c\left\|f^{\prime \prime \prime}\right\|_{\infty} .$$

  comment

• 1.I.6F

  Numerical Analysis

  Part IB, 2005

  Determine the Cholesky factorization (without pivoting) of the matrix

  $$A=\left[\begin{array}{ccc} 2 & -4 & 2 \\ -4 & 10+\lambda & 2+3 \lambda \\ 2 & 2+3 \lambda & 23+9 \lambda \end{array}\right]$$

  where $\lambda$ is a real parameter. Hence, find the range of values of $\lambda$ for which the matrix $A$ is positive definite.

  comment

• 2.II.18F

  Numerical Analysis

  Part IB, 2005

  (a) Let $\left\{Q_{n}\right\}_{n \geqslant 0}$ be a set of polynomials orthogonal with respect to some inner product $(\cdot, \cdot)$ in the interval $[a, b]$. Write explicitly the least-squares approximation to $f \in C[a, b]$ by an $n$ th-degree polynomial in terms of the polynomials $\left\{Q_{n}\right\}_{n \geqslant 0}$.

  (b) Let an inner product be defined by the formula

  $$(g, h)=\int_{-1}^{1}\left(1-x^{2}\right)^{-\frac{1}{2}} g(x) h(x) d x$$

  Determine the $n$th degree polynomial approximation of $f(x)=\left(1-x^{2}\right)^{\frac{1}{2}}$ with respect to this inner product as a linear combination of the underlying orthogonal polynomials.

  comment

• 3.II.19F

  Numerical Analysis

  Part IB, 2005

  Given real $\mu \neq 0$, we consider the matrix

  $$A=\left[\begin{array}{cccc} \frac{1}{\mu} & 1 & 0 & 0 \\ -1 & \frac{1}{\mu} & 1 & 0 \\ 0 & -1 & \frac{1}{\mu} & 1 \\ 0 & 0 & -1 & \frac{1}{\mu} \end{array}\right]$$

  Construct the Jacobi and Gauss-Seidel iteration matrices originating in the solution of the linear system $A x=b$.

  Determine the range of real $\mu \neq 0$ for which each iterative procedure converges.

  comment

• 4.I.8F

  Numerical Analysis

  Part IB, 2005

  Define Gaussian quadrature.

  Evaluate the coefficients of the Gaussian quadrature of the integral

  $$\int_{-1}^{1}\left(1-x^{2}\right) f(x) d x$$

  which uses two function evaluations.

  comment

• 1.I.6D

  Numerical Analysis

  Part IB, 2006

  (a) Perform the LU-factorization with column pivoting of the matrix

  $$A=\left[\begin{array}{rrr} 2 & 1 & 1 \\ 4 & 1 & 0 \\ -2 & 2 & 1 \end{array}\right]$$

  (b) Explain briefly why every nonsingular matrix $A$ admits an LU-factorization with column pivoting.

  comment

• 2.II.18D

  Numerical Analysis

  Part IB, 2006

  (a) For a positive weight function $w$, let

  $$\int_{-1}^{1} f(x) w(x) d x \approx \sum_{i=0}^{n} a_{i} f\left(x_{i}\right)$$

  be the corresponding Gaussian quadrature with $n+1$ nodes. Prove that all the coefficients $a_{i}$ are positive.

  (b) The integral

  $$I(f)=\int_{-1}^{1} f(x) w(x) d x$$

  is approximated by a quadrature

  $$I_{n}(f)=\sum_{i=0}^{n} a_{i}^{(n)} f\left(x_{i}^{(n)}\right)$$

  which is exact on polynomials of degree $\leqslant n$ and has positive coefficients $a_{i}^{(n)}$. Prove that, for any $f$ continuous on $[-1,1]$, the quadrature converges to the integral, i.e.,

  $$\left|I(f)-I_{n}(f)\right| \rightarrow 0 \quad \text { as } \quad n \rightarrow \infty$$

  [You may use the Weierstrass theorem: for any $f$ continuous on $[-1,1]$, and for any $\epsilon>0$, there exists a polynomial $Q$ of degree $n=n(\epsilon, f)$ such that $\left.\max _{x \in[-1,1]}|f(x)-Q(x)|<\epsilon .\right]$

  comment

• 3.II.19D

  Numerical Analysis

  Part IB, 2006

  (a) Define the QR factorization of a rectangular matrix and explain how it can be used to solve the least squares problem of finding an $x^{*} \in \mathbb{R}^{n}$ such that

  $$\left\|A x^{*}-b\right\|=\min _{x \in \mathbb{R}^{n}}\|A x-b\|, \quad \text { where } \quad A \in \mathbb{R}^{m \times n}, \quad b \in \mathbb{R}^{m}, \quad m \geqslant n,$$

  and the norm is the Euclidean distance $\|y\|=\sqrt{\sum_{i=1}^{m}\left|y_{i}\right|^{2}}$.

  (b) Define a Householder transformation (reflection) $H$ and prove that $H$ is an orthogonal matrix.

  (c) Using Householder reflection, solve the least squares problem for the case

  $$A=\left[\begin{array}{rr} 2 & 4 \\ 1 & -1 \\ 2 & 1 \end{array}\right], \quad b=\left[\begin{array}{l} 1 \\ 5 \\ 1 \end{array}\right]$$

  and find the value of $\left\|A x^{*}-b\right\|=\min _{x \in \mathbb{R}^{2}}\|A x-b\|$.

  comment

• 4.I.8D

  Numerical Analysis

  Part IB, 2006

  (a) Given the data

  \begin{tabular}{c|r|r|r|r} $x_{i}$ & $-1$ & 0 & 1 & 3 \ \hline$f\left(x_{i}\right)$ & $-7$ & $-3$ & $-3$ & 9 \end{tabular}

  find the interpolating cubic polynomial $p \in \mathcal{P}_{3}$ in the Newton form, and transform it to the power form.

  (b) We add to the data one more value $f\left(x_{i}\right)$ at $x_{i}=2$. Find the power form of the interpolating quartic polynomial $q \in \mathcal{P}_{4}$ to the extended data

  \begin{tabular}{c|c|r|r|r|r} $x_{i}$ & $-1$ & 0 & 1 & 2 & 3 \ \hline$f\left(x_{i}\right)$ & $-7$ & $-3$ & $-3$ & $-7$ & 9 \end{tabular}$.$

  comment

• 1.I.6F

  Numerical Analysis

  Part IB, 2007

  Solve the least squares problem

  $$\left[\begin{array}{ll} 1 & 3 \\ 0 & 2 \\ 0 & 2 \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} x_{1} \\ x_{2} \end{array}\right]=\left[\begin{array}{r} 4 \\ 1 \\ 4 \\ -1 \end{array}\right]$$

  using $Q R$ method with Householder transformation. (A solution using normal equations is not acceptable.)

  comment

• 2.II.18F

  Numerical Analysis

  Part IB, 2007

  For a symmetric, positive definite matrix $A$ with the spectral radius $\rho(A)$, the linear system $A x=b$ is solved by the iterative procedure

  $$x^{(k+1)}=x^{(k)}-\tau\left(A x^{(k)}-b\right), \quad k \geq 0$$

  where $\tau$ is a real parameter. Find the range of $\tau$ that guarantees convergence of $x^{(k)}$ to the exact solution for any choice of $x^{(0)}$.

  comment

• 3.II.19F

  Numerical Analysis

  Part IB, 2007

  Prove that the monic polynomials $Q_{n}, n \geq 0$, orthogonal with respect to a given weight function $w(x)>0$ on $[a, b]$, satisfy the three-term recurrence relation

  $$Q_{n+1}(x)=\left(x-a_{n}\right) Q_{n}(x)-b_{n} Q_{n-1}(x), \quad n \geq 0$$

  where $Q_{-1}(x) \equiv 0, Q_{0}(x) \equiv 1$. Express the values $a_{n}$ and $b_{n}$ in terms of $Q_{n}$ and $Q_{n-1}$ and show that $b_{n}>0$.

  comment

• 4.I.8F

  Numerical Analysis

  Part IB, 2007

  Given $f \in C^{3}[0,2]$, we approximate $f^{\prime}(0)$ by the linear combination

  $$\mu(f)=-\frac{3}{2} f(0)+2 f(1)-\frac{1}{2} f(2)$$

  Using the Peano kernel theorem, determine the least constant $c$ in the inequality

  $$\left|f^{\prime}(0)-\mu(f)\right| \leq c\left\|f^{\prime \prime \prime}\right\|_{\infty},$$

  and give an example of $f$ for which the inequality turns into equality.

  comment

• 1.I.6D

  Numerical Analysis

  Part IB, 2008

  Show that if $A=L D L^{T}$, where $L \in \mathbb{R}^{m \times m}$ is a lower triangular matrix with all elements on the main diagonal being unity and $D \in \mathbb{R}^{m \times m}$ is a diagonal matrix with positive elements, then $A$ is positive definite. Find $L$ and the corresponding $D$ when

  $$A=\left[\begin{array}{rrr} 1 & -1 & 2 \\ -1 & 3 & 1 \\ 2 & 1 & 3 \end{array}\right]$$

  comment

• 2.II.18D

  Numerical Analysis

  Part IB, 2008

  (a) A Householder transformation (reflection) is given by

  $$H=I-\frac{2 u u^{T}}{\|u\|^{2}},$$

  where $H \in \mathbb{R}^{m \times m}, u \in \mathbb{R}^{m}$, and $I$ is the $m \times m$ unit matrix and $u$ is a non-zero vector which has norm $\|u\|=\left(\sum_{i=1}^{m} u_{i}^{2}\right)^{1 / 2}$. Show that $H$ is orthogonal.

  (b) Suppose that $A \in \mathbb{R}^{m \times n}, x \in \mathbb{R}^{n}$ and $b \in \mathbb{R}^{m}$ with $n<m$. Show that if $x$ minimises $\|A x-b\|^{2}$ then it also minimises $\|Q A x-Q b\|^{2}$, where $Q$ is an arbitrary $m \times m$ orthogonal matrix.

  (c) Using Householder reflection, find the $x$ that minimises $\|A x-b\|^{2}$ when

  $$A=\left[\begin{array}{ll} 1 & 2 \\ 0 & 4 \\ 0 & 2 \\ 0 & 4 \end{array}\right] \quad b=\left[\begin{array}{r} 1 \\ 1 \\ 2 \\ -1 \end{array}\right]$$

  comment

• 3.II.19D

  Numerical Analysis

  Part IB, 2008

  Starting from the Taylor formula for $f(x) \in C^{k+1}[a, b]$ with an integral remainder term, show that the error of an approximant $L(f)$ can be written in the form (Peano kernel theorem)

  $$L(f)=\frac{1}{k !} \int_{a}^{b} K(\theta) f^{(k+1)}(\theta) d \theta,$$

  when $L(f)$, which is identically zero if $f(x)$ is a polynomial of degree $k$, satisfies conditions that you should specify. Give an expression for $K(\theta)$.

  Hence determine the minimum value of $c$ in the inequality

  $$|L(f)| \leq c\left\|f^{\prime \prime \prime}\right\|_{\infty}$$

  when

  $$L(f)=f^{\prime}(1)-\frac{1}{2}(f(2)-f(0)) \text { for } f(x) \in C^{3}[0,2]$$

  comment

• 4.I.8D

  Numerical Analysis

  Part IB, 2008

  Show that the Chebyshev polynomials, $T_{n}(x)=\cos \left(n \cos ^{-1} x\right), n=0,1,2, \ldots$ obey the orthogonality relation

  $$\int_{-1}^{1} \frac{T_{n}(x) T_{m}(x)}{\sqrt{1-x^{2}}} d x=\frac{\pi}{2} \delta_{n, m}\left(1+\delta_{n, 0}\right)$$

  State briefly how an optimal choice of the parameters $a_{k}, x_{k}, k=1,2 \ldots n$ is made in the Gaussian quadrature formula

  $$\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^{2}}} d x \sim \sum_{k=1}^{n} a_{k} f\left(x_{k}\right)$$

  Find these parameters for the case $n=3$.

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2009

  The real non-singular matrix $A \in \mathbb{R}^{m \times m}$ is written in the form $A=A_{D}+A_{U}+A_{L}$, where the matrices $A_{D}, A_{U}, A_{L} \in \mathbb{R}^{m \times m}$ are diagonal and non-singular, strictly uppertriangular and strictly lower-triangular respectively.

  Given $b \in \mathbb{R}^{m}$, the Jacobi iteration for solving $A x=b$ is

  $$A_{D} x_{n}=-\left(A_{U}+A_{L}\right) x_{n-1}+b, \quad n=1,2 \ldots$$

  where the $n$th iterate is $x_{n} \in \mathbb{R}^{m}$. Show that the iteration converges to the solution $x$ of $A x=b$, independent of the starting choice $x_{0}$, if and only if the spectral radius $\rho(H)$ of the matrix $H=-A_{D}^{-1}\left(A_{U}+A_{L}\right)$ is less than 1 .

  Hence find the range of values of the real number $\mu$ for which the iteration will converge when

  $$A=\left[\begin{array}{rrr} 1 & 0 & -\mu \\ -\mu & 3 & -\mu \\ -4 \mu & 0 & 4 \end{array}\right]$$

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2009

  Suppose that $w(x)>0$ for all $x \in(a, b)$. The weights $b_{1}, \ldots, b_{n}$ and nodes $x_{1}, \ldots, x_{n}$ are chosen so that the Gaussian quadrature formula

  $$\int_{a}^{b} w(x) f(x) d x \sim \sum_{k=1}^{n} b_{k} f\left(x_{k}\right)$$

  is exact for every polynomial of degree $2 n-1$. Show that the $b_{i}, i=1, \ldots, n$ are all positive.

  When $w(x)=1+x^{2}, a=-1$ and $b=1$, the first three underlying orthogonal polynomials are $p_{0}(x)=1, p_{1}(x)=x$, and $p_{2}(x)=x^{2}-2 / 5$. Find $x_{1}, x_{2}$ and $b_{1}, b_{2}$ when $n=2$.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2009

  The real orthogonal matrix $\Omega^{[p, q]} \in \mathbb{R}^{m \times m}$ with $1 \leqslant p<q \leqslant m$ is a Givens rotation with rotation angle $\theta$. Write down the form of $\Omega^{[p, q]}$.

  Show that for any matrix $A \in \mathbb{R}^{m \times m}$ it is possible to choose $\theta$ such that the matrix $\Omega^{[p, q]} A$ satisfies $\left(\Omega^{[p, q]} A\right)_{q, j}=0$ for any $j$, where $1 \leqslant j \leqslant m$.

  Let

  $$A=\left[\begin{array}{ccc} 1 & 3 & 2 \\ 1 & 4 & 4 \\ \sqrt{2} & 7 / \sqrt{2} & 4 \sqrt{2} \end{array}\right]$$

  By applying a sequence of Givens rotations of the form $\Omega^{[1,3]} \Omega^{[1,2]}$, chosen to reduce the elements in the first column below the main diagonal to zero, find a factorisation of the matrix $A \in \mathbb{R}^{3 \times 3}$ of the form $A=Q R$, where $Q \in \mathbb{R}^{3 \times 3}$ is an orthogonal matrix and $R \in \mathbb{R}^{3 \times 3}$ is an upper-triangular matrix for which the leading non-zero element in each row is positive.

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2009

  Starting from Taylor's theorem with integral form of the remainder, prove the Peano kernel theorem: the error of an approximant $L(f)$ applied to $f(x) \in C^{k+1}[a, b]$ can be written in the form

  $$L(f)=\frac{1}{k !} \int_{a}^{b} K(\theta) f^{(k+1)}(\theta) d \theta$$

  You should specify the form of $K(\theta)$. Here it is assumed that $L(f)$ is identically zero when $f(x)$ is a polynomial of degree $k$. State any other necessary conditions.

  Setting $a=0$ and $b=2$, find $K(\theta)$ and show that it is negative for $0<\theta<2$ when

  $$L(f)=\int_{0}^{2} f(x) d x-\frac{1}{3}(f(0)+4 f(1)+f(2)) \quad \text { for } \quad f(x) \in C^{4}[0,2] .$$

  Hence determine the minimum value of $\rho$ for which

  $$|L(f)| \leqslant \rho\left\|f^{(4)}\right\|_{\infty}$$

  holds for all $f(x) \in C^{4}[0,2]$.

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2010

  Obtain the Cholesky decompositions of

  $$H_{3}=\left(\begin{array}{ccc} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{array}\right), \quad H_{4}=\left(\begin{array}{cccc} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \lambda \end{array}\right) .$$

  What is the minimum value of $\lambda$ for $H_{4}$ to be positive definite? Verify that if $\lambda=\frac{1}{7}$ then $H_{4}$ is positive definite.

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2010

  Suppose $x_{0}, x_{1}, \ldots, x_{n} \in[a, b] \subset \mathbf{R}$ are pointwise distinct and $f(x)$ is continuous on $[a, b]$. For $k=1,2, \ldots, n$ define

  $$I_{0, k}(x)=\frac{f\left(x_{0}\right)\left(x_{k}-x\right)-f\left(x_{k}\right)\left(x_{0}-x\right)}{x_{k}-x_{0}}$$

  and for $k=2,3, \ldots, n$

  $$I_{0,1, \ldots, k-2, k-1, k}(x)=\frac{I_{0,1, \ldots, k-2, k-1}(x)\left(x_{k}-x\right)-I_{0,1, \ldots, k-2, k}(x)\left(x_{k-1}-x\right)}{x_{k}-x_{k-1}}$$

  Show that $I_{0,1}, \ldots, k-2, k-1, k(x)$ is a polynomial of order $k$ which interpolates $f(x)$ at $x_{0}, x_{1}, \ldots, x_{k}$.

  Given $x_{k}=\{-1,0,2,5\}$ and $f\left(x_{k}\right)=\{33,5,9,1335\}$, determine the interpolating polynomial.

  comment

• Paper 1, Section II, 18C

  Numerical Analysis

  Part IB, 2010

  Let

  $$\langle f, g\rangle=\int_{-\infty}^{\infty} e^{-x^{2}} f(x) g(x) d x$$

  be an inner product. The Hermite polynomials $H_{n}(x), n=0,1,2, \ldots$ are polynomials in $x$ of degree $n$ with leading term $2^{n} x^{n}$ which are orthogonal with respect to the inner product, with

  $$\left\langle H_{m}, H_{n}\right\rangle= \begin{cases}\gamma_{m}>0 & \text { if } m=n \\ 0 & \text { otherwise }\end{cases}$$

  and $H_{0}(x)=1$. Find a three-term recurrence relation which is satisfied by $H_{n}(x)$ and $\gamma_{n}$ for $n=1,2,3$. [You may assume without proof that

  $$\left.\langle 1,1\rangle=\sqrt{\pi}, \quad\langle x, x\rangle=\frac{1}{2} \sqrt{\pi}, \quad\left\langle x^{2}, x^{2}\right\rangle=\frac{3}{4} \sqrt{\pi}, \quad\left\langle x^{3}, x^{3}\right\rangle=\frac{15}{8} \sqrt{\pi} .\right]$$

  Next let $x_{0}, x_{1}, \ldots, x_{k}$ be the $k+1$ distinct zeros of $H_{k+1}(x)$ and for $i, j=0,1, \ldots, k$ define the Lagrangian polynomials

  $$L_{i}(x)=\prod_{j \neq i} \frac{x-x_{j}}{x_{i}-x_{j}}$$

  associated with these points. Prove that $\left\langle L_{i}, L_{j}\right\rangle=0$ if $i \neq j$.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2010

  Consider the initial value problem for an autonomous differential equation

  $$y^{\prime}(t)=f(y(t)), \quad y(0)=y_{0} \text { given }$$

  and its approximation on a grid of points $t_{n}=n h, n=0,1,2, \ldots$. Writing $y_{n}=y\left(t_{n}\right)$, it is proposed to use one of two Runge-Kutta schemes defined by

  $$y_{n+1}=y_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right)$$

  where $k_{1}=h f\left(y_{n}\right)$ and

  $$k_{2}= \begin{cases}h f\left(y_{n}+k_{1}\right) & \text { scheme I } \\ h f\left(y_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right)\right) & \text { scheme II }\end{cases}$$

  What is the order of each scheme? Determine the $A$-stability of each scheme.

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2010

  Define the QR factorization of an $m \times n$ matrix $A$ and explain how it can be used to solve the least squares problem of finding the $x^{*} \in \mathbf{R}^{n}$ which minimises $\|A x-b\|$ where $b \in \mathbf{R}^{m}, m>n$, and the norm is the Euclidean one.

  Define a Householder (reflection) transformation $H$ and show that it is an orthogonal matrix.

  Using a Householder reflection, solve the least squares problem for

  $$A=\left(\begin{array}{rrr} 2 & 4 & 7 \\ 0 & 3 & -1 \\ 0 & 0 & 2 \\ 0 & 0 & 1 \\ 0 & 0 & -2 \end{array}\right), \quad b=\left(\begin{array}{r} 9 \\ -7 \\ 3 \\ 1 \\ -1 \end{array}\right)$$

  giving both $x^{*}$ and $\left\|A x^{*}-b\right\|$.

  comment

• Paper 1, Section I, B

  Numerical Analysis

  Part IB, 2011

  Orthogonal monic polynomials $p_{0}, p_{1}, \ldots, p_{n}, \ldots$ are defined with respect to the inner product $\langle p, q\rangle=\int_{-1}^{1} w(x) p(x) q(x) d x$, where $p_{n}$ is of degree $n$. Show that such polynomials obey a three-term recurrence relation

  $$p_{n+1}(x)=\left(x-\alpha_{n}\right) p_{n}(x)-\beta_{n} p_{n-1}(x)$$

  for appropriate choices of $\alpha_{n}$ and $\beta_{n}$.

  Now suppose that $w(x)$ is an even function of $x$. Show that the $p_{n}$ are even or odd functions of $x$ according to whether $n$ is even or odd.

  comment

• Paper 4, Section I, B

  Numerical Analysis

  Part IB, 2011

  Consider the multistep method for numerical solution of the differential equation $\mathbf{y}^{\prime}=\mathbf{f}(t, \mathbf{y})$ :

  $$\sum_{l=0}^{s} \rho_{l} \mathbf{y}_{n+l}=h \sum_{l=0}^{s} \sigma_{l} \mathbf{f}\left(t_{n+l}, \mathbf{y}_{n+l}\right), \quad n=0,1, \ldots$$

  What does it mean to say that the method is of order $p$, and that the method is convergent?

  Show that the method is of order $p$ if

  $$\sum_{l=0}^{s} \rho_{l}=0, \quad \sum_{l=0}^{s} l^{k} \rho_{l}=k \sum_{l=0}^{s} l^{k-1} \sigma_{l}, \quad k=1,2, \ldots, p$$

  and give the conditions on $\rho(w)=\sum_{l=0}^{s} \rho_{l} w^{l}$ that ensure convergence.

  Hence determine for what values of $\theta$ and the $\sigma_{i}$ the two-step method

  $$\mathbf{y}_{n+2}-(1-\theta) \mathbf{y}_{n+1}-\theta \mathbf{y}_{n}=h\left[\sigma_{0} \mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right)+\sigma_{1} \mathbf{f}\left(t_{n+1}, \mathbf{y}_{n+1}\right)+\sigma_{2} \mathbf{f}\left(t_{n+2}, \mathbf{y}_{n+2}\right)\right]$$

  is (a) convergent, and (b) of order 3 .

  comment

• Paper 1, Section II, B

  Numerical Analysis

  Part IB, 2011

  Consider a function $f(x)$ defined on the domain $x \in[0,1]$. Find constants $\alpha, \beta$, so that for any fixed $\xi \in[0,1]$,

  $$f^{\prime \prime}(\xi)=\alpha f(0)+\beta f^{\prime}(0)+\gamma f(1)$$

  is exactly satisfied for polynomials of degree less than or equal to two.

  By using the Peano kernel theorem, or otherwise, show that

  $$\begin{aligned} f^{\prime}(\xi)-f^{\prime}(0)-\xi(\alpha f(0)&\left.+\beta f^{\prime}(0)+\gamma f(1)\right)=\int_{0}^{\xi}(\xi-\theta) H_{1}(\theta) f^{\prime \prime \prime}(\theta) d \theta \\ &+\int_{0}^{\xi} \theta H_{2}(\theta) f^{\prime \prime \prime}(\theta) d \theta+\int_{\xi}^{1} \xi H_{2}(\theta) f^{\prime \prime \prime}(\theta) d \theta \end{aligned}$$

  where $H_{1}(\theta)=1-(1-\theta)^{2} \geqslant 0, H_{2}(\theta)=-(1-\theta)^{2} \leqslant 0$. Thus show that

  $$\left|f^{\prime}(\xi)-f^{\prime}(0)-\xi\left(\alpha f(0)+\beta f^{\prime}(0)+\gamma f(1)\right)\right| \leqslant\left.\frac{1}{6}\left(2 \xi-3 \xi^{2}+4 \xi^{3}-\xi^{4}\right)|| f^{\prime \prime \prime}\right|_{\infty} .$$

  comment

• Paper 2, Section II, B

  Numerical Analysis

  Part IB, 2011

  What is the $Q R$-decomposition of a matrix A? Explain how to construct the matrices $Q$ and $R$ by the Gram-Schmidt procedure, and show how the decomposition can be used to solve the matrix equation $A \mathbf{x}=\mathbf{b}$ when $A$ is a square matrix.

  Why is this procedure not useful for numerical decomposition of large matrices? Give a brief description of an alternative procedure using Givens rotations.

  Find a $Q R$-decomposition for the matrix

  $$\mathrm{A}=\left[\begin{array}{rrrr} 3 & 4 & 7 & 13 \\ -6 & -8 & -8 & -12 \\ 3 & 4 & 7 & 11 \\ 0 & 2 & 5 & 7 \end{array}\right]$$

  Is your decomposition unique? Use the decomposition you have found to solve the equation

  $$A x=\left[\begin{array}{c} 4 \\ 6 \\ 2 \\ 9 \end{array}\right]$$

  comment

• Paper 3, Section II, B

  Numerical Analysis

  Part IB, 2011

  A Gaussian quadrature formula provides an approximation to the integral

  $$\int_{-1}^{1}\left(1-x^{2}\right) f(x) d x \approx \sum_{k=1}^{\nu} b_{k} f\left(c_{k}\right)$$

  which is exact for all $f(x)$ that are polynomials of degree $\leqslant(2 \nu-1)$.

  Write down explicit expressions for the $b_{k}$ in terms of integrals, and explain why it is necessary that the $c_{k}$ are the zeroes of a (monic) polynomial $p_{\nu}$ of degree $\nu$ that satisfies $\int_{-1}^{1}\left(1-x^{2}\right) p_{\nu}(x) q(x) d x=0$ for any polynomial $q(x)$ of degree less than $\nu .$

  The first such polynomials are $p_{0}=1, p_{1}=x, p_{2}=x^{2}-1 / 5, p_{3}=x^{3}-3 x / 7$. Show that the Gaussian quadrature formulae for $\nu=2,3$ are

  $$\begin{array}{ll} \nu=2: & \frac{2}{3}\left[f\left(-\frac{1}{\sqrt{5}}\right)+f\left(\frac{1}{\sqrt{5}}\right)\right] \\ \nu=3: & \frac{14}{45}\left[f\left(-\sqrt{\frac{3}{7}}\right)+f\left(\sqrt{\frac{3}{7}}\right)\right]+\frac{32}{45} f(0) \end{array}$$

  Verify the result for $\nu=3$ by considering $f(x)=1, x^{2}, x^{4}$.

  comment

• Paper 1, Section I, D

  Numerical Analysis

  Part IB, 2012

  $$A=\left[\begin{array}{cccc} 1 & a & a^{2} & a^{3} \\ a^{3} & 1 & a & a^{2} \\ a^{2} & a^{3} & 1 & a \\ a & a^{2} & a^{3} & 1 \end{array}\right], \quad b=\left[\begin{array}{l} \gamma \\ 0 \\ 0 \\ 0 \end{array}\right], \quad \gamma=1-a^{4} \neq 0$$

  Find the LU factorization of the matrix $A$ and use it to solve the system $A x=b$ via forward and backward substitution. [Other methods of solution are not acceptable.]

  comment

• Paper 4, Section I, D

  Numerical Analysis

  Part IB, 2012

  State the Dahlquist equivalence theorem regarding convergence of a multistep method.

  The multistep method, with a real parameter $a$,

  $$y_{n+3}+(2 a-3)\left(y_{n+2}-y_{n+1}\right)-y_{n}=h a\left(f_{n+2}-f_{n+1}\right)$$

  is of order 2 for any $a$, and also of order 3 for $a=6$. Determine all values of $a$ for which the method is convergent, and find the order of convergence.

  comment

• Paper 1, Section II, D

  Numerical Analysis

  Part IB, 2012

  For a numerical method for solving $y^{\prime}=f(t, y)$, define the linear stability domain, and state when such a method is A-stable.

  Determine all values of the real parameter $a$ for which the Runge-Kutta method

  $$\begin{aligned} k_{1} &=f\left(t_{n}+\left(\frac{1}{2}-a\right) h, y_{n}+\left(\frac{1}{4} h k_{1}+\left(\frac{1}{4}-a\right) h k_{2}\right)\right) \\ k_{2} &=f\left(t_{n}+\left(\frac{1}{2}+a\right) h, y_{n}+\left(\left(\frac{1}{4}+a\right) h k_{1}+\frac{1}{4} h k_{2}\right)\right) \\ y_{n+1} &=y_{n}+\frac{1}{2} h\left(k_{1}+k_{2}\right) \end{aligned}$$

  is A-stable.

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part IB, 2012

  Define the QR factorization of an $m \times n$ matrix $A$ and explain how it can be used to solve the least squares problem of finding the vector $x^{*} \in \mathbb{R}^{n}$ which minimises $\left\|A x^{*}-b\right\|$, where $b \in \mathbb{R}^{m}, m>n$, and the norm is the Euclidean one.

  Define a Householder transformation $H$ and show that it is an orthogonal matrix.

  Using a Householder transformation, solve the least squares problem for

  $$A=\left[\begin{array}{rrr} 1 & -1 & 5 \\ 0 & 1 & 5 \\ 0 & 0 & 3 \\ 0 & 0 & 4 \end{array}\right], \quad b=\left[\begin{array}{r} 1 \\ 2 \\ -1 \\ 2 \end{array}\right]$$

  giving both $x^{*}$ and $\left\|A x^{*}-b\right\|$.

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part IB, 2012

  Let $\left\{P_{n}\right\}_{n=0}^{\infty}$ be the sequence of monic polynomials of degree $n$ orthogonal on the interval $[-1,1]$ with respect to the weight function $w$.

  Prove that each $P_{n}$ has $n$ distinct zeros in the interval $(-1,1)$.

  Let $P_{0}(x)=1, P_{1}(x)=x-a_{1}$, and let $P_{n}$ satisfy the following three-term recurrence relation:

  $$P_{n}(x)=\left(x-a_{n}\right) P_{n-1}(x)-b_{n}^{2} P_{n-2}(x), \quad n \geqslant 2$$

  Set

  $$A_{n}=\left[\begin{array}{ccccc} a_{1} & b_{2} & 0 & \cdots & 0 \\ b_{2} & a_{2} & b_{3} & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & b_{n-1} & a_{n-1} & b_{n} \\ 0 & \cdots & 0 & b_{n} & a_{n} \end{array}\right]$$

  Prove that $P_{n}(x)=\operatorname{det}\left(x I-A_{n}\right), n \geqslant 1$, and deduce that all the eigenvalues of $A_{n}$ are distinct and reside in $(-1,1)$.

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2013

  Determine the nodes $x_{1}, x_{2}$ of the two-point Gaussian quadrature

  $$\int_{0}^{1} f(x) w(x) d x \approx a_{1} f\left(x_{1}\right)+a_{2} f\left(x_{2}\right), \quad w(x)=x$$

  and express the coefficients $a_{1}, a_{2}$ in terms of $x_{1}, x_{2}$. [You don't need to find numerical values of the coefficients.]

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2013

  For a continuous function $f$, and $k+1$ distinct points $\left\{x_{0}, x_{1}, \ldots, x_{k}\right\}$, define the divided difference $f\left[x_{0}, \ldots, x_{k}\right]$ of order $k$.

  Given $n+1$ points $\left\{x_{0}, x_{1}, \ldots, x_{n}\right\}$, let $p_{n} \in \mathbb{P}_{n}$ be the polynomial of degree $n$ that interpolates $f$ at these points. Prove that $p_{n}$ can be written in the Newton form

  $$p_{n}(x)=f\left(x_{0}\right)+\sum_{k=1}^{n} f\left[x_{0}, \ldots, x_{k}\right] \prod_{i=0}^{k-1}\left(x-x_{i}\right)$$

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part IB, 2013

  Define the QR factorization of an $m \times n$ matrix $A$ and explain how it can be used to solve the least squares problem of finding the vector $x^{*} \in \mathbb{R}^{n}$ which minimises $\|A x-b\|$, where $b \in \mathbb{R}^{m}, m>n$, and the norm is the Euclidean one.

  Define a Givens rotation $\Omega^{[p, q]}$ and show that it is an orthogonal matrix.

  Using a Givens rotation, solve the least squares problem for

  $$A=\left[\begin{array}{lll} 2 & 1 & 1 \\ 0 & 4 & 1 \\ 0 & 3 & 2 \\ 0 & 0 & 0 \end{array}\right], \quad b=\left[\begin{array}{l} 2 \\ 3 \\ 1 \\ 2 \end{array}\right]$$

  giving both $x^{*}$ and $\left\|A x^{*}-b\right\|$.

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2013

  $$f^{\prime}(0) \approx a_{0} f(0)+a_{1} f(1)+a_{2} f(2)=: \lambda(f)$$

  be a formula of numerical differentiation which is exact on polynomials of degree 2 , and let

  $$e(f)=f^{\prime}(0)-\lambda(f)$$

  be its error.

  Find the values of the coefficients $a_{0}, a_{1}, a_{2}$.

  Using the Peano kernel theorem, find the least constant $c$ such that, for all functions $f \in C^{3}[0,2]$, we have

  $$|e(f)| \leqslant c\left\|f^{\prime \prime \prime}\right\|_{\infty} .$$

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2013

  Explain briefly what is meant by the convergence of a numerical method for solving the ordinary differential equation

  $$y^{\prime}(t)=f(t, y), \quad t \in[0, T], \quad y(0)=y_{0} .$$

  Prove from first principles that if the function $f$ is sufficiently smooth and satisfies the Lipschitz condition

  $$|f(t, x)-f(t, y)| \leqslant L|x-y|, \quad x, y \in \mathbb{R}, \quad t \in[0, T],$$

  for some $L>0$, then the backward Euler method

  $$y_{n+1}=y_{n}+h f\left(t_{n+1}, y_{n+1}\right)$$

  converges and find the order of convergence.

  Find the linear stability domain of the backward Euler method.

  comment

• Paper 1, Section I, $\mathbf{6 C}$

  Numerical Analysis

  Part IB, 2014

  (i) A general multistep method for the numerical approximation to the scalar differential equation $y^{\prime}=f(t, y)$ is given by

  $$\sum_{\ell=0}^{s} \rho_{\ell} y_{n+\ell}=h \sum_{\ell=0}^{s} \sigma_{\ell} f_{n+\ell}, \quad n=0,1, \ldots$$

  where $f_{n+\ell}=f\left(t_{n+\ell}, y_{n+\ell}\right)$. Show that this method is of order $p \geqslant 1$ if and only if

  $$\rho\left(\mathrm{e}^{z}\right)-z \sigma\left(\mathrm{e}^{z}\right)=\mathcal{O}\left(z^{p+1}\right) \quad \text { as } \quad z \rightarrow 0$$

  where

  $$\rho(w)=\sum_{\ell=0}^{s} \rho_{\ell} w^{\ell} \quad \text { and } \quad \sigma(w)=\sum_{\ell=0}^{s} \sigma_{\ell} w^{\ell}$$

  (ii) A particular three-step implicit method is given by

  $$y_{n+3}+(a-1) y_{n+1}-a y_{n}=h\left(f_{n+3}+\sum_{\ell=0}^{2} \sigma_{\ell} f_{n+\ell}\right)$$

  where the $\sigma_{\ell}$ are chosen to make the method third order. [The $\sigma_{\ell}$ need not be found.] For what values of $a$ is the method convergent?

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2014

  Consider the quadrature given by

  $$\int_{0}^{\pi} w(x) f(x) d x \approx \sum_{k=1}^{\nu} b_{k} f\left(c_{k}\right)$$

  for $\nu \in \mathbb{N}$, disjoint $c_{k} \in(0, \pi)$ and $w>0$. Show that it is not possible to make this quadrature exact for all polynomials of order $2 \nu$.

  For the case that $\nu=2$ and $w(x)=\sin x$, by considering orthogonal polynomials find suitable $b_{k}$ and $c_{k}$ that make the quadrature exact on cubic polynomials.

  [Hint: $\int_{0}^{\pi} x^{2} \sin x d x=\pi^{2}-4$ and $\int_{0}^{\pi} x^{3} \sin x d x=\pi^{3}-6 \pi .$ ]

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part IB, 2014

  Define a Householder transformation $\mathrm{H}$ and show that it is an orthogonal matrix. Briefly explain how these transformations can be used for QR factorisation of an $m \times n$ matrix.

  Using Householder transformations, find a QR factorisation of

  $$A=\left[\begin{array}{rrr} 2 & 5 & 4 \\ 2 & 5 & 1 \\ -2 & 1 & 5 \\ 2 & -1 & 16 \end{array}\right]$$

  Using this factorisation, find the value of $\lambda$ for which

  $$A x=\left[\begin{array}{c} 1+\lambda \\ 2 \\ 3 \\ 4 \end{array}\right]$$

  has a unique solution $x \in \mathbb{R}^{3}$.

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2014

  A Runge-Kutta scheme is given by

  $$k_{1}=h f\left(y_{n}\right), \quad k_{2}=h f\left(y_{n}+\left[(1-a) k_{1}+a k_{2}\right]\right), \quad y_{n+1}=y_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right)$$

  for the solution of an autonomous differential equation $y^{\prime}=f(y)$, where $a$ is a real parameter. What is the order of the scheme? Identify all values of $a$ for which the scheme is A-stable. Determine the linear stability domain for this range.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2014

  A linear functional acting on $f \in C^{k+1}[a, b]$ is approximated using a linear scheme $L(f)$. The approximation is exact when $f$ is a polynomial of degree $k$. The error is given by $\lambda(f)$. Starting from the Taylor formula for $f(x)$ with an integral remainder term, show that the error can be written in the form

  $$\lambda(f)=\frac{1}{k !} \int_{a}^{b} K(\theta) f^{(k+1)}(\theta) d \theta$$

  subject to a condition on $\lambda$ that you should specify. Give an expression for $K(\theta)$.

  Find $c_{0}, c_{1}$ and $c_{3}$ such that the approximation scheme

  $$f^{\prime \prime}(2) \approx c_{0} f(0)+c_{1} f(1)+c_{3} f(3)$$

  is exact for all $f$ that are polynomials of degree 2 . Assuming $f \in C^{3}[0,3]$, apply the Peano kernel theorem to the error. Find and sketch $K(\theta)$ for $k=2$.

  Find the minimum values for the constants $r$ and $s$ for which

  $$|\lambda(f)| \leqslant r\left\|f^{(3)}\right\|_{1} \quad \text { and } \quad|\lambda(f)| \leqslant s\left\|f^{(3)}\right\|_{\infty}$$

  and show explicitly that both error bounds hold for $f(x)=x^{3}$.

  comment

• Paper 1, Section I, 6D

  Numerical Analysis

  Part IB, 2015

  Let

  $$A=\left[\begin{array}{rrrr} 1 & 4 & 3 & 2 \\ 4 & 17 & 13 & 11 \\ 3 & 13 & 13 & 12 \\ 2 & 11 & 12 & \lambda \end{array}\right], \quad b=\left[\begin{array}{l} 1 \\ 1 \\ 3 \\ 2 \end{array}\right]$$

  where $\lambda$ is a real parameter. Find the $L U$ factorization of the matrix $A$. Give the constraint on $\lambda$ for A to be positive definite.

  For $\lambda=18$, use this factorization to solve the system $A x=b$ via forward and backward substitution.

  comment

• Paper 4, Section I, D

  Numerical Analysis

  Part IB, 2015

  Given $n+1$ distinct points $\left\{x_{0}, x_{1}, \ldots, x_{n}\right\}$, let $p_{n} \in \mathbb{P}_{n}$ be the real polynomial of degree $n$ that interpolates a continuous function $f$ at these points. State the Lagrange interpolation formula.

  Prove that $p_{n}$ can be written in the Newton form

  $$p_{n}(x)=f\left(x_{0}\right)+\sum_{k=1}^{n} f\left[x_{0}, \ldots, x_{k}\right] \prod_{i=0}^{k-1}\left(x-x_{i}\right)$$

  where $f\left[x_{0}, \ldots, x_{k}\right]$ is the divided difference, which you should define. [An explicit expression for the divided difference is not required.]

  Explain why it can be more efficient to use the Newton form rather than the Lagrange formula.

  comment

• Paper 1, Section II, 18D

  Numerical Analysis

  Part IB, 2015

  Determine the real coefficients $b_{1}, b_{2}, b_{3}$ such that

  $$\int_{-2}^{2} f(x) d x=b_{1} f(-1)+b_{2} f(0)+b_{3} f(1)$$

  is exact when $f(x)$ is any real polynomial of degree 2 . Check explicitly that the quadrature is exact for $f(x)=x^{2}$ with these coefficients.

  State the Peano kernel theorem and define the Peano kernel $K(\theta)$. Use this theorem to show that if $f \in C^{3}[-2,2]$, and $b_{1}, b_{2}, b_{3}$ are chosen as above, then

  $$\left|\int_{-2}^{2} f(x) d x-b_{1} f(-1)-b_{2} f(0)-b_{3} f(1)\right| \leqslant \frac{4}{9} \max _{\xi \in[-2,2]}\left|f^{(3)}(\xi)\right|$$

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part IB, 2015

  Define the QR factorization of an $m \times n$ matrix $A$. Explain how it can be used to solve the least squares problem of finding the vector $x^{*} \in \mathbb{R}^{n}$ which minimises $\|\mathrm{A} x-b\|$, where $b \in \mathbb{R}^{m}, m>n$, and $\|\cdot\|$ is the Euclidean norm.

  Explain how to construct $Q$ and $R$ by the Gram-Schmidt procedure. Why is this procedure not useful for numerical factorization of large matrices?

  Let

  $$A=\left[\begin{array}{rrr} 5 & 6 & -14 \\ 5 & 4 & 4 \\ -5 & 2 & -8 \\ 5 & 12 & -18 \end{array}\right], \quad b=\left[\begin{array}{l} 1 \\ 1 \\ 1 \\ 0 \end{array}\right]$$

  Using the Gram-Schmidt procedure find a QR decomposition of A. Hence solve the least squares problem giving both $x^{*}$ and $\left\|\mathrm{A} x^{*}-b\right\|$.

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part IB, 2015

  Define the linear stability domain for a numerical method to solve $y^{\prime}=f(t, y)$. What is meant by an $A$-stable method? Briefly explain the relevance of these concepts in the numerical solution of ordinary differential equations.

  Consider

  $$y_{n+1}=y_{n}+h\left[\theta f\left(t_{n}, y_{n}\right)+(1-\theta) f\left(t_{n+1}, y_{n+1}\right)\right]$$

  where $\theta \in[0,1]$. What is the order of this method?

  Find the linear stability domain of this method. For what values of $\theta$ is the method A-stable?

  comment

• Paper 1, Section I, D

  Numerical Analysis

  Part IB, 2016

  (a) What are real orthogonal polynomials defined with respect to an inner product $\langle\cdot, \cdot\rangle ?$ What does it mean for such polynomials to be monic?

  (b) Real monic orthogonal polynomials, $p_{n}(x)$, of degree $n=0,1,2, \ldots$, are defined with respect to the inner product,

  $$\langle p, q\rangle=\int_{-1}^{1} w(x) p(x) q(x) d x$$

  where $w(x)$ is a positive weight function. Show that such polynomials obey the three-term recurrence relation,

  $$p_{n+1}(x)=\left(x-\alpha_{n}\right) p_{n}(x)-\beta_{n} p_{n-1}(x),$$

  for appropriate $\alpha_{n}$ and $\beta_{n}$ which should be given in terms of inner products.

  comment

• Paper 4, Section I, D

  Numerical Analysis

  Part IB, 2016

  (a) Define the linear stability domain for a numerical method to solve $\mathbf{y}^{\prime}=\mathbf{f}(t, \mathbf{y})$. What is meant by an A-stable method?

  (b) A two-stage Runge-Kutta scheme is given by

  $$\mathbf{k}_{1}=\mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right), \quad \mathbf{k}_{2}=\mathbf{f}\left(t_{n}+\frac{h}{2}, \mathbf{y}_{n}+\frac{h}{2} \mathbf{k}_{1}\right), \quad \mathbf{y}_{n+1}=\mathbf{y}_{n}+h \mathbf{k}_{2},$$

  where $h$ is the step size and $t_{n}=n h$. Show that the order of this scheme is at least two. For this scheme, find the intersection of the linear stability domain with the real axis. Hence show that this method is not A-stable.

  comment

• Paper 1, Section II, D

  Numerical Analysis

  Part IB, 2016

  (a) Consider a method for numerically solving an ordinary differential equation (ODE) for an initial value problem, $\mathbf{y}^{\prime}=\mathbf{f}(t, \mathbf{y})$. What does it mean for a method to converge over $t \in[0, T]$ where $T \in \mathbb{R}$ ? What is the definition of the order of a method?

  (b) A general multistep method for the numerical solution of an ODE is

  $$\sum_{l=0}^{s} \rho_{l} \mathbf{y}_{n+l}=h \sum_{l=0}^{s} \sigma_{l} \mathbf{f}\left(t_{n+l}, \mathbf{y}_{n+l}\right), \quad n=0,1, \ldots$$

  where $s$ is a fixed positive integer. Show that this method is at least of order $p \geqslant 1$ if and only if

  $$\sum_{l=0}^{s} \rho_{l}=0 \quad \text { and } \quad \sum_{l=0}^{s} l^{k} \rho_{l}=k \sum_{l=0}^{s} l^{k-1} \sigma_{l}, \quad k=1, \ldots, p$$

  (c) State the Dahlquist equivalence theorem regarding the convergence of a multistep method.

  (d) Consider the multistep method,

  $$\mathbf{y}_{n+2}+\theta \mathbf{y}_{n+1}+a \mathbf{y}_{n}=h\left[\sigma_{0} \mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right)+\sigma_{1} \mathbf{f}\left(t_{n+1}, \mathbf{y}_{n+1}\right)+\sigma_{2} \mathbf{f}\left(t_{n+2}, \mathbf{y}_{n+2}\right)\right]$$

  Determine the values of $\sigma_{i}$ and $a$ (in terms of the real parameter $\theta$ ) such that the method is at least third order. For what values of $\theta$ does the method converge?

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part IB, 2016

  (a) Determine real quadratic functions $a(x), b(x), c(x)$ such that the interpolation formula,

  $$f(x) \approx a(x) f(0)+b(x) f(2)+c(x) f(3),$$

  is exact when $f(x)$ is any real polynomial of degree 2 .

  (b) Use this formula to construct approximations for $f(5)$ and $f^{\prime}(1)$ which are exact when $f(x)$ is any real polynomial of degree 2 . Calculate these approximations for $f(x)=x^{3}$ and comment on your answers.

  (c) State the Peano kernel theorem and define the Peano kernel $K(\theta)$. Use this theorem to find the minimum values of the constants $\alpha$ and $\beta$ such that

  $$\left|f(1)-\frac{1}{3}[f(0)+3 f(2)-f(3)]\right| \leqslant \alpha \max _{\xi \in[0,3]}\left|f^{(2)}(\xi)\right|$$

  and

  $$\left|f(1)-\frac{1}{3}[f(0)+3 f(2)-f(3)]\right| \leqslant \beta\left\|f^{(2)}\right\|_{1}$$

  where $f \in C^{2}[0,3]$. Check that these inequalities hold for $f(x)=x^{3}$.

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part IB, 2016

  (a) Define a Givens rotation $\Omega^{[p, q]} \in \mathbb{R}^{m \times m}$ and show that it is an orthogonal matrix.

  (b) Define a QR factorization of a matrix $A \in \mathbb{R}^{m \times n}$ with $m \geqslant n$. Explain how Givens rotations can be used to find $Q \in \mathbb{R}^{m \times m}$ and $R \in \mathbb{R}^{m \times n}$.

  (c) Let

  $$\mathbf{A}=\left[\begin{array}{ccc} 3 & 1 & 1 \\ 0 & 4 & 1 \\ 0 & 3 & 2 \\ 0 & 0 & 3 / 4 \end{array}\right], \quad \mathbf{b}=\left[\begin{array}{c} 98 / 25 \\ 25 \\ 25 \\ 0 \end{array}\right]$$

  (i) Find a QR factorization of $A$ using Givens rotations.

  (ii) Hence find the vector $\mathbf{x}^{*} \in \mathbb{R}^{3}$ which minimises $\|A \mathbf{x}-\mathbf{b}\|$, where $\|\cdot\|$ is the Euclidean norm. What is $\left\|\mathrm{A} \mathbf{x}^{*}-\mathbf{b}\right\|$ ?

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2017

  Given $n+1$ real points $x_{0}<x_{1}<\cdots<x_{n}$, define the Lagrange cardinal polynomials $\ell_{i}(x), i=0,1, \ldots, n$. Let $p(x)$ be the polynomial of degree $n$ that interpolates the function $f \in C^{n}\left[x_{0}, x_{n}\right]$ at these points. Express $p(x)$ in terms of the values $f_{i}=f\left(x_{i}\right)$, $i=0,1, \ldots, n$ and the Lagrange cardinal polynomials.

  Define the divided difference $f\left[x_{0}, x_{1}, \ldots, x_{n}\right]$ and give an expression for it in terms of $f_{0}, f_{1}, \ldots, f_{n}$ and $x_{0}, x_{1}, \ldots, x_{n}$. Prove that

  $$f\left[x_{0}, x_{1}, \ldots, x_{n}\right]=\frac{1}{n !} f^{(n)}(\xi)$$

  for some number $\xi \in\left[x_{0}, x_{n}\right]$.

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2017

  For the matrix

  $$A=\left[\begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & 5 & 5 & 5 \\ 1 & 5 & 14 & 14 \\ 1 & 5 & 14 & \lambda \end{array}\right]$$

  find a factorization of the form

  $$A=L D L^{\top} \text {, }$$

  where $D$ is diagonal and $L$ is lower triangular with ones on its diagonal.

  For what values of $\lambda$ is $A$ positive definite?

  In the case $\lambda=30$ find the Cholesky factorization of $A$.

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part IB, 2017

  A three-stage explicit Runge-Kutta method for solving the autonomous ordinary differential equation

  $$\frac{d y}{d t}=f(y)$$

  is given by

  $$y_{n+1}=y_{n}+h\left(b_{1} k_{1}+b_{2} k_{2}+b_{3} k_{3}\right),$$

  where

  $$\begin{aligned} &k_{1}=f\left(y_{n}\right) \\ &k_{2}=f\left(y_{n}+h a_{1} k_{1}\right) \\ &k_{3}=f\left(y_{n}+h\left(a_{2} k_{1}+a_{3} k_{2}\right)\right) \end{aligned}$$

  and $h>0$ is the time-step. Derive sufficient conditions on the coefficients $b_{1}, b_{2}, b_{3}, a_{1}$, $a_{2}$ and $a_{3}$ for the method to be of third order.

  Assuming that these conditions hold, verify that $-\frac{5}{2}$ belongs to the linear stability domain of the method.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2017

  Define the linear least-squares problem for the equation $A \mathbf{x}=\mathbf{b}$, where $A$ is an $m \times n$ matrix with $m>n, \mathbf{b} \in \mathbb{R}^{m}$ is a given vector and $\mathbf{x} \in \mathbb{R}^{n}$ is an unknown vector.

  If $A=Q R$, where $Q$ is an orthogonal matrix and $R$ is an upper triangular matrix in standard form, explain why the least-squares problem is solved by minimizing the Euclidean norm $\left\|R \mathbf{x}-Q^{\top} \mathbf{b}\right\|$.

  Using the method of Householder reflections, find a QR factorization of the matrix

  $$A=\left[\begin{array}{rrr} 1 & 3 & 3 \\ 1 & 3 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & -1 \end{array}\right]$$

  Hence find the solution of the least-squares problem in the case

  $$\mathbf{b}=\left[\begin{array}{r} 1 \\ 1 \\ 3 \\ -1 \end{array}\right]$$

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2017

  Let $p_{n} \in \mathbb{P}_{n}$ be the $n$th monic orthogonal polynomial with respect to the inner product

  $$\langle f, g\rangle=\int_{a}^{b} w(x) f(x) g(x) d x$$

  on $C[a, b]$, where $w$ is a positive weight function.

  Prove that, for $n \geqslant 1, p_{n}$ has $n$ distinct zeros in the interval $(a, b)$.

  Let $c_{1}, c_{2}, \ldots, c_{n} \in[a, b]$ be $n$ distinct points. Show that the quadrature formula

  $$\int_{a}^{b} w(x) f(x) d x \approx \sum_{i=1}^{n} b_{i} f\left(c_{i}\right)$$

  is exact for all $f \in \mathbb{P}_{n-1}$ if the weights $b_{i}$ are chosen to be

  $$b_{i}=\int_{a}^{b} w(x) \prod_{\substack{j=1 \\ j \neq i}}^{n} \frac{x-c_{j}}{c_{i}-c_{j}} d x$$

  Show further that the quadrature formula is exact for all $f \in \mathbb{P}_{2 n-1}$ if the nodes $c_{i}$ are chosen to be the zeros of $p_{n}$ (Gaussian quadrature). [Hint: Write $f$ as $q p_{n}+r$, where $q, r \in \mathbb{P}_{n-1}$.]

  Use the Peano kernel theorem to write an integral expression for the approximation error of Gaussian quadrature for sufficiently differentiable functions. (You should give a formal expression for the Peano kernel but are not required to evaluate it.)

  comment

• Paper 1, Section I, D

  Numerical Analysis

  Part IB, 2018

  The Trapezoidal Rule for solving the differential equation

  $$y^{\prime}(t)=f(t, y), \quad t \in[0, T], \quad y(0)=y_{0}$$

  is defined by

  $$y_{n+1}=y_{n}+\frac{1}{2} h\left[f\left(t_{n}, y_{n}\right)+f\left(t_{n+1}, y_{n+1}\right)\right]$$

  where $h=t_{n+1}-t_{n}$.

  Determine the minimum order of convergence $k$ of this rule for general functions $f$ that are sufficiently differentiable. Show with an explicit example that there is a function $f$ for which the local truncation error is $A h^{k+1}$ for some constant $A$.

  comment

• Paper 4 , Section I, D

  Numerical Analysis

  Part IB, 2018

  $$A=\left[\begin{array}{cccc} 1 & 2 & 1 & 2 \\ 2 & 5 & 5 & 6 \\ 1 & 5 & 13 & 14 \\ 2 & 6 & 14 & \lambda \end{array}\right], \quad b=\left[\begin{array}{l} 1 \\ 3 \\ 7 \\ \mu \end{array}\right]$$

  where $\lambda$ and $\mu$ are real parameters. Find the $L U$ factorisation of the matrix $A$. For what values of $\lambda$ does the equation $A x=b$ have a unique solution for $x$ ?

  For $\lambda=20$, use the $L U$ decomposition with forward and backward substitution to determine a value for $\mu$ for which a solution to $A x=b$ exists. Find the most general solution to the equation in this case.

  comment

• Paper 1, Section II, D

  Numerical Analysis

  Part IB, 2018

  Show that if $\mathbf{u} \in \mathbb{R}^{m} \backslash\{\mathbf{0}\}$ then the $m \times m$ matrix transformation

  $$H_{\mathbf{u}}=I-2 \frac{\mathbf{u} \mathbf{u}^{\top}}{\|\mathbf{u}\|^{2}}$$

  is orthogonal. Show further that, for any two vectors $\mathbf{a}, \mathbf{b} \in \mathbb{R}^{m}$ of equal length,

  $$H_{\mathbf{a}-\mathbf{b}} \mathbf{a}=\mathbf{b} .$$

  Explain how to use such transformations to convert an $m \times n$ matrix $A$ with $m \geqslant n$ into the form $A=Q R$, where $Q$ is an orthogonal matrix and $R$ is an upper-triangular matrix, and illustrate the method using the matrix

  $$A=\left[\begin{array}{rrr} 1 & -1 & 4 \\ 1 & 4 & -2 \\ 1 & 4 & 2 \\ 1 & -1 & 0 \end{array}\right]$$

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part IB, 2018

  Taylor's theorem for functions $f \in C^{k+1}[a, b]$ is given in the form

  $$f(x)=f(a)+(x-a) f^{\prime}(a)+\cdots+\frac{(x-a)^{k}}{k !} f^{(k)}(a)+R(x) .$$

  Use integration by parts to show that

  $$R(x)=\frac{1}{k !} \int_{a}^{x}(x-\theta)^{k} f^{(k+1)}(\theta) d \theta$$

  Let $\lambda_{k}$ be a linear functional on $C^{k+1}[a, b]$ such that $\lambda_{k}[p]=0$ for $p \in \mathbb{P}_{k}$. Show that

  $$\lambda_{k}[f]=\frac{1}{k !} \int_{a}^{b} K(\theta) f^{(k+1)}(\theta) d \theta$$

  where the Peano kernel function $K(\theta)=\lambda_{k}\left[(x-\theta)_{+}^{k}\right] . \quad[$ You may assume that the functional commutes with integration over a fixed interval.]

  The error in the mid-point rule for numerical quadrature on $[0,1]$ is given by

  $$e[f]=\int_{0}^{1} f(x) d x-f\left(\frac{1}{2}\right)$$

  Show that $e[p]=0$ if $p$ is a linear polynomial. Find the Peano kernel function corresponding to $e$ explicitly and verify the formula ( $\dagger$ ) in the case $f(x)=x^{2}$.

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part IB, 2018

  Show that the recurrence relation

  $$\begin{aligned} p_{0}(x) &=1 \\ p_{n+1}(x) &=q_{n+1}(x)-\sum_{k=0}^{n} \frac{\left\langle q_{n+1}, p_{k}\right\rangle}{\left\langle p_{k}, p_{k}\right\rangle} p_{k}(x) \end{aligned}$$

  where $\langle\cdot, \cdot\rangle$ is an inner product on real polynomials, produces a sequence of orthogonal, monic, real polynomials $p_{n}(x)$ of degree exactly $n$ of the real variable $x$, provided that $q_{n}$ is a monic, real polynomial of degree exactly $n$.

  Show that the choice $q_{n+1}(x)=x p_{n}(x)$ leads to a three-term recurrence relation of the form

  $$\begin{aligned} p_{0}(x) &=1 \\ p_{1}(x) &=x-\alpha_{0} \\ p_{n+1}(x) &=\left(x-\alpha_{n}\right) p_{n}(x)-\beta_{n} p_{n-1}(x) \end{aligned}$$

  where $\alpha_{n}$ and $\beta_{n}$ are constants that should be determined in terms of the inner products $\left\langle p_{n}, p_{n}\right\rangle,\left\langle p_{n-1}, p_{n-1}\right\rangle$ and $\left\langle p_{n}, x p_{n}\right\rangle$.

  Use this recurrence relation to find the first four monic Legendre polynomials, which correspond to the inner product defined by

  $$\langle p, q\rangle \equiv \int_{-1}^{1} p(x) q(x) d x$$

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2019

  Let $[a, b]$ be the smallest interval that contains the $n+1$ distinct real numbers $x_{0}, x_{1}, \ldots, x_{n}$, and let $f$ be a continuous function on that interval.

  Define the divided difference $f\left[x_{0}, x_{1}, \ldots, x_{m}\right]$ of degree $m \leqslant n$.

  Prove that the polynomial of degree $n$ that interpolates the function $f$ at the points $x_{0}, x_{1}, \ldots, x_{n}$ is equal to the Newton polynomial

  $$p_{n}(x)=f\left[x_{0}\right]+f\left[x_{0}, x_{1}\right]\left(x-x_{0}\right)+\cdots+f\left[x_{0}, x_{1}, \ldots, x_{n}\right] \prod_{i=0}^{n-1}\left(x-x_{i}\right)$$

  Prove the recursive formula

  $$f\left[x_{0}, x_{1}, \ldots, x_{m}\right]=\frac{f\left[x_{1}, x_{2}, \ldots, x_{m}\right]-f\left[x_{0}, x_{1}, \ldots, x_{m-1}\right]}{x_{m}-x_{0}}$$

  for $1 \leqslant m \leqslant n .$

  comment

• Paper 4, Section I, C

  Numerical Analysis

  Part IB, 2019

  Calculate the $L U$ factorization of the matrix

  $$A=\left(\begin{array}{rrrr} 3 & 2 & -3 & -3 \\ 6 & 3 & -7 & -8 \\ 3 & 1 & -6 & -4 \\ -6 & -3 & 9 & 6 \end{array}\right)$$

  Use this to evaluate $\operatorname{det}(A)$ and to solve the equation

  $$A \mathbf{x}=\mathbf{b}$$

  with

  $$\mathbf{b}=\left(\begin{array}{r} 3 \\ 3 \\ -1 \\ -3 \end{array}\right)$$

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part IB, 2019

  (a) An $s$-step method for solving the ordinary differential equation

  $$\frac{d \mathbf{y}}{d t}=\mathbf{f}(t, \mathbf{y})$$

  is given by

  $$\sum_{l=0}^{s} \rho_{l} \mathbf{y}_{n+l}=h \sum_{l=0}^{s} \sigma_{l} \mathbf{f}\left(t_{n+l}, \mathbf{y}_{n+l}\right), \quad n=0,1, \ldots$$

  where $\rho_{l}$ and $\sigma_{l}(l=0,1, \ldots, s)$ are constant coefficients, with $\rho_{s}=1$, and $h$ is the time-step. Prove that the method is of order $p \geqslant 1$ if and only if

  $$\rho\left(e^{z}\right)-z \sigma\left(e^{z}\right)=O\left(z^{p+1}\right)$$

  as $z \rightarrow 0$, where

  $$\rho(w)=\sum_{l=0}^{s} \rho_{l} w^{l}, \quad \sigma(w)=\sum_{l=0}^{s} \sigma_{l} w^{l}$$

  (b) Show that the Adams-Moulton method

  $$\mathbf{y}_{n+2}=\mathbf{y}_{n+1}+\frac{h}{12}\left(5 \mathbf{f}\left(t_{n+2}, \mathbf{y}_{n+2}\right)+8 \mathbf{f}\left(t_{n+1}, \mathbf{y}_{n+1}\right)-\mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right)\right)$$

  is of third order and convergent.

  [You may assume the Dahlquist equivalence theorem if you state it clearly.]

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part IB, 2019

  (a) Let $w(x)$ be a positive weight function on the interval $[a, b]$. Show that

  $$\langle f, g\rangle=\int_{a}^{b} f(x) g(x) w(x) d x$$

  defines an inner product on $C[a, b]$.

  (b) Consider the sequence of polynomials $p_{n}(x)$ defined by the three-term recurrence relation

  $$p_{n+1}(x)=\left(x-\alpha_{n}\right) p_{n}(x)-\beta_{n} p_{n-1}(x), \quad n=1,2, \ldots$$

  where

  $$p_{0}(x)=1, \quad p_{1}(x)=x-\alpha_{0},$$

  and the coefficients $\alpha_{n}$ (for $\left.n \geqslant 0\right)$ and $\beta_{n}$ (for $\left.n \geqslant 1\right)$ are given by

  $$\alpha_{n}=\frac{\left\langle p_{n}, x p_{n}\right\rangle}{\left\langle p_{n}, p_{n}\right\rangle}, \quad \beta_{n}=\frac{\left\langle p_{n}, p_{n}\right\rangle}{\left\langle p_{n-1}, p_{n-1}\right\rangle}$$

  Prove that this defines a sequence of monic orthogonal polynomials on $[a, b]$.

  (c) The Hermite polynomials $H e_{n}(x)$ are orthogonal on the interval $(-\infty, \infty)$ with weight function $e^{-x^{2} / 2}$. Given that

  $$H e_{n}(x)=(-1)^{n} e^{x^{2} / 2} \frac{d^{n}}{d x^{n}}\left(e^{-x^{2} / 2}\right)$$

  deduce that the Hermite polynomials satisfy a relation of the form $(*)$ with $\alpha_{n}=0$ and $\beta_{n}=n$. Show that $\left\langle H e_{n}, H e_{n}\right\rangle=n ! \sqrt{2 \pi}$.

  (d) State, without proof, how the properties of the Hermite polynomial $\operatorname{He}_{N}(x)$, for some positive integer $N$, can be used to estimate the integral

  $$\int_{-\infty}^{\infty} f(x) e^{-x^{2} / 2} d x$$

  where $f(x)$ is a given function, by the method of Gaussian quadrature. For which polynomials is the quadrature formula exact?

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2019

  Define the linear least squares problem for the equation

  $$A \mathbf{x}=\mathbf{b}$$

  where $A$ is a given $m \times n$ matrix with $m>n, \mathbf{b} \in \mathbb{R}^{m}$ is a given vector and $\mathbf{x} \in \mathbb{R}^{n}$ is an unknown vector.

  Explain how the linear least squares problem can be solved by obtaining a $Q R$ factorization of the matrix $A$, where $Q$ is an orthogonal $m \times m$ matrix and $R$ is an uppertriangular $m \times n$ matrix in standard form.

  Use the Gram-Schmidt method to obtain a $Q R$ factorization of the matrix

  $$A=\left(\begin{array}{lll} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 0 & 0 \end{array}\right)$$

  and use it to solve the linear least squares problem in the case

  $$\mathbf{b}=\left(\begin{array}{l} 1 \\ 2 \\ 3 \\ 6 \end{array}\right)$$

  comment

• Paper 1, Section I, C

  Numerical Analysis

  Part IB, 2020

  (a) Find an $L U$ factorisation of the matrix

  $$A=\left[\begin{array}{cccc} 1 & 1 & 0 & 3 \\ 0 & 2 & 2 & 12 \\ 0 & 5 & 7 & 32 \\ 3 & -1 & -1 & -10 \end{array}\right]$$

  where the diagonal elements of $L$ are $L_{11}=L_{44}=1, L_{22}=L_{33}=2$.

  (b) Use this factorisation to solve the linear system $A \mathbf{x}=\mathbf{b}$, where

  $$\mathbf{b}=\left[\begin{array}{c} -3 \\ -12 \\ -30 \\ 13 \end{array}\right]$$

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part IB, 2020

  (a) Given a set of $n+1$ distinct real points $x_{0}, x_{1}, \ldots, x_{n}$ and real numbers $f_{0}, f_{1}, \ldots, f_{n}$, show that the interpolating polynomial $p_{n} \in \mathbb{P}_{n}[x], p_{n}\left(x_{i}\right)=f_{i}$, can be written in the form

  $$p_{n}(x)=\sum_{k=0}^{n} a_{k} \prod_{j=0, j \neq k}^{n} \frac{x-x_{j}}{x_{k}-x_{j}}, \quad x \in \mathbb{R}$$

  where the coefficients $a_{k}$ are to be determined.

  (b) Consider the approximation of the integral of a function $f \in C[a, b]$ by a finite sum,

  $$\int_{a}^{b} f(x) d x \approx \sum_{k=0}^{s-1} w_{k} f\left(c_{k}\right)$$

  where the weights $w_{0}, \ldots, w_{s-1}$ and nodes $c_{0}, \ldots, c_{s-1} \in[a, b]$ are independent of $f$. Derive the expressions for the weights $w_{k}$ that make the approximation ( 1$)$ exact for $f$ being any polynomial of degree $s-1$, i.e. $f \in \mathbb{P}_{s-1}[x]$.

  Show that by choosing $c_{0}, \ldots, c_{s-1}$ to be zeros of the polynomial $q_{s}(x)$ of degree $s$, one of a sequence of orthogonal polynomials defined with respect to the scalar product

  $$\langle u, v\rangle=\int_{a}^{b} u(x) v(x) d x$$

  the approximation (1) becomes exact for $f \in \mathbb{P}_{2 s-1}[x]$ (i.e. for all polynomials of degree $2 s-1)$.

  (c) On the interval $[a, b]=[-1,1]$ the scalar product (2) generates orthogonal polynomials given by

  $$q_{n}(x)=\frac{1}{2^{n} n !} \frac{d^{n}}{d x^{n}}\left(x^{2}-1\right)^{n}, \quad n=0,1,2, \ldots$$

  Find the values of the nodes $c_{k}$ for which the approximation (1) is exact for all polynomials of degree 7 (i.e. $f \in \mathbb{P}_{7}[x]$ ) but no higher.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part IB, 2020

  Consider a multistep method for numerical solution of the differential equation $\mathbf{y}^{\prime}=\mathbf{f}(t, \mathbf{y})$ :

  $$\mathbf{y}_{n+2}-\mathbf{y}_{n+1}=h\left[(1+\alpha) \mathbf{f}\left(t_{n+2}, \mathbf{y}_{n+2}\right)+\beta \mathbf{f}\left(t_{n+1}, \mathbf{y}_{n+1}\right)-(\alpha+\beta) \mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right)\right],$$

  where $n=0,1, \ldots$, and $\alpha$ and $\beta$ are constants.

  (a) Define the order of a method for numerically solving an ODE.

  (b) Show that in general an explicit method of the form $(*)$ has order 1 . Determine the values of $\alpha$ and $\beta$ for which this multistep method is of order 3 .

  (c) Show that the multistep method (*) is convergent.

  comment

• Paper 1, Section I, B

  Numerical Analysis

  Part IB, 2021

  Prove, from first principles, that there is an algorithm that can determine whether any real symmetric matrix $A \in \mathbb{R}^{n \times n}$ is positive definite or not, with the computational cost (number of arithmetic operations) bounded by $\mathcal{O}\left(n^{3}\right)$.

  [Hint: Consider the LDL decomposition.]

  comment

• Paper 4, Section I, B

  Numerical Analysis

  Part IB, 2021

  (a) Given the data $f(0)=0, f(1)=4, f(2)=2, f(3)=8$, find the interpolating cubic polynomial $p_{3} \in \mathbb{P}_{3}[x]$ in the Newton form.

  (b) We add to the data one more value, $f(-2)=10$. Find the interpolating quartic polynomial $p_{4} \in \mathbb{P}_{4}[x]$ for the extended data in the Newton form.

  comment

• Paper 1, Section II, B

  Numerical Analysis

  Part IB, 2021

  For the ordinary differential equation

  $$\boldsymbol{y}^{\prime}=\boldsymbol{f}(t, \boldsymbol{y}), \quad \boldsymbol{y}(0)=\tilde{\boldsymbol{y}}_{0}, \quad t \geqslant 0$$

  where $\boldsymbol{y}(t) \in \mathbb{R}^{N}$ and the function $\boldsymbol{f}: \mathbb{R} \times \mathbb{R}^{N} \rightarrow \mathbb{R}^{N}$ is analytic, consider an explicit one-step method described as the mapping

  $$\boldsymbol{y}_{n+1}=\boldsymbol{y}_{n}+h \varphi\left(t_{n}, \boldsymbol{y}_{n}, h\right)$$

  Here $\varphi: \mathbb{R}_{+} \times \mathbb{R}^{N} \times \mathbb{R}_{+} \rightarrow \mathbb{R}^{N}, n=0,1, \ldots$ and $t_{n}=n h$ with time step $h>0$, producing numerical approximations $\boldsymbol{y}_{n}$ to the exact solution $\boldsymbol{y}\left(t_{n}\right)$ of equation $(*)$, with $\boldsymbol{y}_{0}$ being the initial value of the numerical solution.

  (i) Define the local error of a one-step method.

  (ii) Let $\|\cdot\|$ be a norm on $\mathbb{R}^{N}$ and suppose that

  $$\|\boldsymbol{\varphi}(t, \boldsymbol{u}, h)-\boldsymbol{\varphi}(t, \boldsymbol{v}, h)\| \leqslant L\|\boldsymbol{u}-\boldsymbol{v}\|,$$

  for all $h>0, t \in \mathbb{R}, \boldsymbol{u}, \boldsymbol{v} \in \mathbb{R}^{N}$, where $L$ is some positive constant. Let $t^{*}>0$ be given and $\boldsymbol{e}_{0}=\boldsymbol{y}_{0}-\boldsymbol{y}(0)$ denote the initial error (potentially non-zero). Show that if the local error of the one-step method ( $\uparrow$ ) is $\mathcal{O}\left(h^{p+1}\right)$, then

  $$\max _{n=0, \ldots,\left\lfloor t^{*} / h\right\rfloor}\left\|\boldsymbol{y}_{n}-\boldsymbol{y}(n h)\right\| \leqslant e^{t^{*} L}\left\|\boldsymbol{e}_{0}\right\|+\mathcal{O}\left(h^{p}\right), \quad h \rightarrow 0$$

  (iii) Let $N=1$ and consider equation $(*)$ where $f$ is time-independent satisfying $|f(u)-f(v)| \leqslant K|u-v|$ for all $u, v \in \mathbb{R}$, where $K$ is a positive constant. Consider the one-step method given by

  $$y_{n+1}=y_{n}+\frac{1}{4} h\left(k_{1}+3 k_{2}\right), \quad k_{1}=f\left(y_{n}\right), \quad k_{2}=f\left(y_{n}+\frac{2}{3} h k_{1}\right) .$$

  Use part (ii) to show that for this method we have that equation (††) holds (with a potentially different constant $L$ ) for $p=2$.

  comment

• Paper 2, Section II, 17B

  Numerical Analysis

  Part IB, 2021

  (a) Define Householder reflections and show that a real Householder reflection is symmetric and orthogonal. Moreover, show that if $H, A \in \mathbb{R}^{n \times n}$, where $H$ is a Householder reflection and $A$ is a full matrix, then the computational cost (number of arithmetic operations) of computing $H A H^{-1}$ can be $\mathcal{O}\left(n^{2}\right)$ operations, as opposed to $\mathcal{O}\left(n^{3}\right)$ for standard matrix products.

  (b) Show that for any $A \in \mathbb{R}^{n \times n}$ there exists an orthogonal matrix $Q \in \mathbb{R}^{n \times n}$ such that

  $$Q A Q^{T}=T=\left[\begin{array}{ccccc} t_{1,1} & t_{1,2} & t_{1,3} & \cdots & t_{1, n} \\ t_{2,1} & t_{2,2} & t_{2,3} & \cdots & t_{2, n} \\ 0 & t_{3,2} & t_{3,3} & \cdots & t_{3, n} \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & t_{n, n-1} & t_{n, n} \end{array}\right]$$

  In particular, $T$ has zero entries below the first subdiagonal. Show that one can compute such a $T$ and $Q$ (they may not be unique) using $\mathcal{O}\left(n^{3}\right)$ arithmetic operations.

  [Hint: Multiply A from the left and right with Householder reflections.]

  comment

• Paper 3, Section II, B

  Numerical Analysis

  Part IB, 2021

  The functions $p_{0}, p_{1}, p_{2}, \ldots$ are generated by the formula

  $$p_{n}(x)=(-1)^{n} x^{-1 / 2} e^{x} \frac{d^{n}}{d x^{n}}\left(x^{n+1 / 2} e^{-x}\right), \quad 0 \leqslant x<\infty$$

  (a) Show that $p_{n}(x)$ is a monic polynomial of degree $n$. Write down the explicit forms of $p_{0}(x), p_{1}(x), p_{2}(x)$.

  (b) Demonstrate the orthogonality of these polynomials with respect to the scalar product

  $$\langle f, g\rangle=\int_{0}^{\infty} x^{1 / 2} e^{-x} f(x) g(x) d x$$

  i.e. that $\left\langle p_{n}, p_{m}\right\rangle=0$ for $m \neq n$, and show that

  $$\left\langle p_{n}, p_{n}\right\rangle=n ! \Gamma\left(n+\frac{3}{2}\right)$$

  where $\Gamma(y)=\int_{0}^{\infty} x^{y-1} e^{-x} d x$.

  (c) Assuming that a three-term recurrence relation in the form

  $$p_{n+1}(x)=\left(x-\alpha_{n}\right) p_{n}(x)-\beta_{n} p_{n-1}(x), \quad n=1,2, \ldots$$

  holds, find the explicit expressions for $\alpha_{n}$ and $\beta_{n}$ as functions of $n$.

  [Hint: you may use the fact that $\Gamma(y+1)=y \Gamma(y) .]$

  comment

• A1.20 B1.20

  Numerical Analysis

  Part II, 2001

  (i) Let $A$ be a symmetric $n \times n$ matrix such that

  $$A_{k, k}>\sum_{\substack{l=1 \\ l \neq k}}^{n}\left|A_{k, l}\right| \quad 1 \leqslant k \leqslant n .$$

  Prove that it is positive definite.

  (ii) Prove that both Jacobi and Gauss-Seidel methods for the solution of the linear system $A \mathrm{x}=\mathbf{b}$, where the matrix $A$ obeys the conditions of (i), converge.

  [You may quote the Householder-John theorem without proof.]

  comment

• A2.19 B2.19

  Numerical Analysis

  Part II, 2001

  (i) Define $m$-step BDF (backward differential formula) methods for the numerical solution of ordinary differential equations and derive explicitly their coefficients.

  (ii) Prove that the linear stability domain of the two-step BDF method includes the interval $(-\infty, 0)$.

  comment

• A3.19 B3.20

  Numerical Analysis

  Part II, 2001

  (i) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}$$

  is discretized by the finite-difference method

  $$u_{m}^{n+1}-\frac{1}{2}(\mu-\alpha)\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right)=u_{m}^{n}+\frac{1}{2}(\mu+\alpha)\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)$$

  where $u_{m}^{n} \approx u(m \Delta x, n \Delta t), \mu=\Delta t /(\Delta x)^{2}$ and $\alpha$ is a constant. Derive the order of magnitude (as a power of $\Delta x$ ) of the local error for different choices of $\alpha$.

  (ii) Investigate the stability of the above finite-difference method for different values of $\alpha$ by the Fourier technique.

  comment

• A4.22 B4.20

  Numerical Analysis

  Part II, 2001

  Write an essay on the computation of eigenvalues and eigenvectors of matrices.

  comment

• A1.20 B1.20

  Numerical Analysis

  Part II, 2002

  (i) Let $A$ be an $n \times n$ symmetric real matrix with distinct eigenvalues $\lambda_{1}, \lambda_{2}, \ldots, \lambda_{n}$ and corresponding eigenvectors $\mathbf{v}_{1}, \mathbf{v}_{2}, \ldots, \mathbf{v}_{n}$, where $\left\|\mathbf{v}_{l}\right\|=1$. Given $\mathbf{x}^{(0)} \in \mathbb{R}^{n},\left\|\mathbf{x}^{(0)}\right\|=1$, the sequence $\mathbf{x}^{(k)}$ is generated in the following manner. We set

  $$\begin{gathered} \mu=\mathbf{x}^{(k) T} A \mathbf{x}^{(k)} \\ \mathbf{y}=(A-\mu I)^{-1} \mathbf{x}^{(k)} \\ \mathbf{x}^{(k+1)}=\frac{\mathbf{y}}{\|\mathbf{y}\|} \end{gathered}$$

  Show that if

  $$\mathbf{x}^{(k)}=c^{-1}\left(\mathbf{v}_{1}+\alpha \sum_{l=2}^{n} d_{l} \mathbf{v}_{l}\right)$$

  where $\alpha$ is a real scalar and $c$ is chosen so that $\left\|\mathbf{x}^{(k)}\right\|=1$, then

  $$\mu=c^{-2}\left(\lambda_{1}+\alpha^{2} \sum_{j=2}^{n} \lambda_{j} d_{j}^{2}\right)$$

  Give an explicit expression for $c$.

  (ii) Use the above result to prove that, if $|\alpha|$ is small,

  $$\mathbf{x}^{(k+1)}=\tilde{c}^{-1}\left(\mathbf{v}_{1}+\alpha^{3} \sum_{l=2}^{n} \tilde{d}_{l} \mathbf{v}_{l}\right)+O\left(\alpha^{4}\right)$$

  and obtain the numbers $\tilde{c}$ and $\tilde{d}_{2}, \ldots, \tilde{d}_{n}$.

  comment

• A2.19 B2.19

  Numerical Analysis

  Part II, 2002

  (i)

  Given the finite-difference method

  $$\sum_{k=-r}^{s} \alpha_{k} u_{m+k}^{n+1}=\sum_{k=-r}^{s} \beta_{k} u_{m+k}^{n}, \quad m, n \in \mathbb{Z}, n \geqslant 0$$

  define

  $$H(z)=\frac{\sum_{k=-r}^{s} \beta_{k} z^{k}}{\sum_{k=-r}^{s} \alpha_{k} z^{k}}$$

  Prove that this method is stable if and only if

  $$\left|H\left(e^{i \theta}\right)\right| \leqslant 1, \quad-\pi \leqslant \theta \leqslant \pi .$$

  [You may quote without proof known properties of the Fourier transform.]

  (ii) Find the range of the parameter $\mu$ such that the method

  $$(1-2 \mu) u_{m-1}^{n+1}+4 \mu u_{m}^{n+1}+(1-2 \mu) u_{m+1}^{n+1}=u_{m-1}^{n}+u_{m+1}^{n}$$

  is stable. Supposing that this method is used to solve the diffusion equation for $u(x, t)$, determine the order of magnitude of the local error as a power of $\Delta x$.

  comment

• A3.19 B3.20

  Numerical Analysis

  Part II, 2002

  (i) Determine the order of the multistep method

  $$\mathbf{y}_{n+2}-(1+\alpha) \mathbf{y}_{n+1}+\alpha \mathbf{y}_{n}=h\left[\frac{1}{12}(5+\alpha) \mathbf{f}_{n+2}+\frac{2}{3}(1-\alpha) \mathbf{f}_{n+1}-\frac{1}{12}(1+5 \alpha) \mathbf{f}_{n}\right]$$

  for the solution of ordinary differential equations for different choices of $\alpha$ in the range $-1 \leqslant \alpha \leqslant 1$.

  (ii) Prove that no such choice of $\alpha$ results in a method whose linear stability domain includes the interval $(-\infty, 0)$.

  comment

• A4.22 B4.20

  Numerical Analysis

  Part II, 2002

  Write an essay on the method of conjugate gradients. You should describe the algorithm, present an analysis of its properties and discuss its advantages.

  [Any theorems quoted should be stated precisely but need not be proved.]

  comment

• A1.20 B1.20

  Numerical Analysis

  Part II, 2003

  (i) The linear algebraic equations $A \mathbf{u}=\mathbf{b}$, where $A$ is symmetric and positive-definite, are solved with the Gauss-Seidel method. Prove that the iteration always converges.

  (ii) The Poisson equation $\nabla^{2} u=f$ is given in the bounded, simply connected domain $\Omega \subseteq \mathbb{R}^{2}$, with zero Dirichlet boundary conditions on $\partial \Omega$. It is approximated by the fivepoint formula

  $$U_{m-1, n}+U_{m, n-1}+U_{m+1, n}+U_{m, n+1}-4 U_{m, n}=(\Delta x)^{2} f_{m, n}$$

  where $U_{m, n} \approx u(m \Delta x, n \Delta x), \quad f_{m, n}=f(m \Delta x, n \Delta x)$, and $(m \Delta x, n \Delta x)$ is in the interior of $\Omega$.

  Assume for the sake of simplicity that the intersection of $\partial \Omega$ with the grid consists only of grid points, so that no special arrangements are required near the boundary. Prove that the method can be written in a vector notation, $A \mathbf{u}=\mathbf{b}$ with a negative-definite matrix $A$.

  comment

• A2.19 B2.19

  Numerical Analysis

  Part II, 2003

  (i) Explain briefly what is meant by the convergence of a numerical method for ordinary differential equations.

  (ii) Suppose the sufficiently-smooth function $f: \mathbb{R} \times \mathbb{R}^{d} \rightarrow \mathbb{R}^{d}$ obeys the Lipschitz condition: there exists $\lambda>0$ such that

  $$\|\mathbf{f}(t, \mathbf{x})-\mathbf{f}(t, \mathbf{y})\| \leqslant \lambda\|\mathbf{x}-\mathbf{y}\|, \quad \mathbf{x}, \mathbf{y} \in \mathbb{R}^{d}, t \geqslant 0$$

  Prove from first principles, without using the Dahlquist equivalence theorem, that the trapezoidal rule

  $$\mathbf{y}_{n+1}=\mathbf{y}_{n}+\frac{1}{2} h\left[\mathbf{f}\left(t_{n}, \mathbf{y}_{n}\right)+\mathbf{f}\left(t_{n+1}, \mathbf{y}_{n+1}\right)\right]$$

  for the solution of the ordinary differential equation

  $$\mathbf{y}^{\prime}=\mathbf{f}(t, \mathbf{y}), \quad t \geqslant 0, \quad \mathbf{y}(0)=\mathbf{y}_{0},$$

  converges.

  comment

• A3.19 B3.20

  Numerical Analysis

  Part II, 2003

  (i) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial}{\partial x}\left(a(x) \frac{\partial u}{\partial x}\right), \quad 0 \leqslant x \leqslant 1, \quad t \geqslant 0$$

  with the initial condition $u(x, 0)=\phi(x), 0 \leqslant x \leqslant 1$ and zero boundary conditions at $x=0$ and $x=1$, is solved by the finite-difference method

  $$\begin{gathered} u_{m}^{n+1}=u_{m}^{n}+\mu\left[a_{m-\frac{1}{2}} u_{m-1}^{n}-\left(a_{m-\frac{1}{2}}+a_{m+\frac{1}{2}}\right) u_{m}^{n}+a_{m+\frac{1}{2}} u_{m+1}^{n}\right] \\ m=1,2, \ldots, N \end{gathered}$$

  where $\mu=\Delta t /(\Delta x)^{2}, \quad \Delta x=\frac{1}{N+1}$ and $u_{m}^{n} \approx u(m \Delta x, n \Delta t), a_{\alpha}=a(\alpha \Delta x)$.

  Assuming sufficient smoothness of the function $a$, and that $\mu$ remains constant as $\Delta x>0$ and $\Delta t>0$ become small, prove that the exact solution satisfies the numerical scheme with error $O\left((\Delta x)^{3}\right)$.

  (ii) For the problem defined in Part (i), assume that there exist $0<a_{-}<a_{+}<\infty$ such that $a_{-} \leqslant a(x) \leqslant a_{+}, \quad 0 \leqslant x \leqslant 1$. Prove that the method is stable for $0<\mu \leqslant 1 /\left(2 a_{+}\right)$.

  [Hint: You may use without proof the Gerschgorin theorem: All the eigenvalues of the matrix $A=\left(a_{k, l}\right)_{k, l=1, \ldots, M}$ are contained in $\bigcup_{k=1}^{m} \mathbb{S}_{k}$, where

  $$\left.\mathbb{S}_{k}=\left\{z \in \mathbb{C}:\left|z-a_{k, k}\right| \leqslant \sum_{\substack{l=1 \\ l \neq k}}^{m}\left|a_{k, l}\right|\right\}, \quad k=1,2, \ldots, m . \quad\right]$$

  comment

• A4.23 B4.20

  Numerical Analysis

  Part II, 2003

  Write an essay on the conjugate gradient method. Your essay should include:

  (a) a statement of the method and a sketch of its derivation;

  (b) discussion, without detailed proofs, but with precise statements of relevant theorems, of the conjugacy of the search directions;

  (c) a description of the standard form of the algorithm;

  (d) discussion of the connection of the method with Krylov subspaces.

  comment

• A1.20 B1.20

  Numerical Analysis

  Part II, 2004

  (i) Define the Backward Difference Formula (BDF) method for ordinary differential equations and derive its two-step version.

  (ii) Prove that the interval $(-\infty, 0)$ belongs to the linear stability domain of the twostep BDF method.

  comment

• A2.19 B2.20

  Numerical Analysis

  Part II, 2004

  (i) The five-point equations, which are obtained when the Poisson equation $\nabla^{2} u=f$ (with Dirichlet boundary conditions) is discretized in a square, are

  $$-u_{m-1, n}-u_{m, n-1}-u_{m+1, n}-u_{m, n+1}+4 u_{m, n}=f_{m, n}, \quad m, n=1,2, \ldots, M,$$

  where $u_{0, n}, u_{M+1, n}, u_{m, 0}, u_{m, M+1}=0$ for all $m, n=1,2, \ldots, M$.

  Formulate the Gauss-Seidel method for the above linear system and prove its convergence. In the proof you should carefully state any theorems you use. [You may use Part (ii) of this question.]

  (ii) By arranging the two-dimensional arrays $\left\{u_{m, n}\right\}_{m, n=1, \ldots, M}$ and $\left\{b_{m, n}\right\}_{m, n=1, \ldots, M}$ into the column vectors $\mathbf{u} \in \mathbb{R}^{M^{2}}$ and $\mathbf{b} \in \mathbb{R}^{M^{2}}$ respectively, the linear system described in Part (i) takes the matrix form $A \mathbf{u}=\mathbf{b}$. Prove that, regardless of the ordering of the points on the grid, the matrix $A$ is symmetric and positive definite.

  comment

• A3.19 B3.20

  Numerical Analysis

  Part II, 2004

  (i) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0 \leqslant x \leqslant 1, \quad t \geqslant 0$$

  with the initial condition $u(x, 0)=\phi(x), 0 \leqslant x \leqslant 1$, and with zero boundary conditions at $x=0$ and $x=1$, can be solved by the method

  $$u_{m}^{n+1}=u_{m}^{n}+\mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right), \quad m=1,2, \ldots, M, \quad n \geqslant 0,$$

  where $\Delta x=1 /(M+1), \mu=\Delta t /(\Delta x)^{2}$, and $u_{m}^{n} \approx u(m \Delta x, n \Delta t)$. Prove that $\mu \leqslant \frac{1}{2}$ implies convergence.

  (ii) By discretizing the same equation and employing the same notation as in Part (i), determine conditions on $\mu>0$ such that the method

  $$\begin{aligned} &\left(\frac{1}{12}-\frac{1}{2} \mu\right) u_{m-1}^{n+1}+\left(\frac{5}{6}+\mu\right) u_{m}^{n+1}+\left(\frac{1}{12}-\frac{1}{2} \mu\right) u_{m+1}^{n+1}= \\ &\left(\frac{1}{12}+\frac{1}{2} \mu\right) u_{m-1}^{n}+\left(\frac{5}{6}-\mu\right) u_{m}^{n}+\left(\frac{1}{12}+\frac{1}{2} \mu\right) u_{m+1}^{n} \end{aligned}$$

  is stable.

  comment

• A4.22 B4.20

  Numerical Analysis

  Part II, 2004

  Write an essay on the method of conjugate gradients. You should define the method, list its main properties and sketch the relevant proof. You should also prove that (in exact arithmetic) the method terminates in a finite number of steps, briefly mention the connection with Krylov subspaces, and describe the approach of preconditioned conjugate gradients.

  comment

• 1.II.38A

  Numerical Analysis

  Part II, 2005

  Let

  $$\frac{\mu}{4} u_{m-1}^{n+1}+u_{m}^{n+1}-\frac{\mu}{4} u_{m+1}^{n+1}=-\frac{\mu}{4} u_{m-1}^{n}+u_{m}^{n}+\frac{\mu}{4} u_{m+1}^{n},$$

  where $n$ is a positive integer and $m$ ranges over all integers, be a finite-difference method for the advection equation

  $$\frac{\partial u}{\partial t}=\frac{\partial u}{\partial x}, \quad-\infty<x<\infty, \quad t \geqslant 0$$

  Here $\mu=\frac{\Delta t}{\Delta x}$ is the Courant number.

  (a) Show that the local error of the method is $O\left((\Delta x)^{3}\right)$.

  (b) Determine the range of $\mu>0$ for which the method is stable.

  comment

• 2.II.38A

  Numerical Analysis

  Part II, 2005

  Define a Krylov subspace $\mathcal{K}_{n}(A, v)$.

  Let $d_{n}$ be the dimension of $\mathcal{K}_{n}(A, v)$. Prove that the sequence $\left\{d_{m}\right\}_{m=1,2, . .}$ increases monotonically. Show that, moreover, there exists an integer $k$ with the following property: $d_{m}=m$ for $m=1,2, \ldots, k$, while $d_{m}=k$ for $m \geqslant k$. Assuming that $A$ has a full set of eigenvectors, show that $k$ is equal to the number of eigenvectors of $A$ required to represent the vector $v$.

  comment

• 3.II.38A

  Numerical Analysis

  Part II, 2005

  Consider the Runge-Kutta method

  $$\begin{aligned} k_{1} &=f\left(y_{n}\right) \\ k_{2} &=f\left(y_{n}+(1-a) h k_{1}+a h k_{2}\right), \\ y_{n+1} &=y_{n}+\frac{h}{2}\left(k_{1}+k_{2}\right) \end{aligned}$$

  for the solution of the scalar ordinary differential equation $y^{\prime}=f(y)$. Here $a$ is a real parameter.

  (a) Determine the order of the method.

  (b) Find the range of values of $a$ for which the method is A-stable.

  comment

• 4.II.39A

  Numerical Analysis

  Part II, 2005

  An $n \times n$ skew-symmetric matrix $A$ is converted into an upper-Hessenberg form $B$, say, by Householder reflections.

  (a) Describe each step of the procedure and observe that $B$ is tridiagonal. Your algorithm should take advantage of the special form of $A$ to reduce the number of computations.

  (b) Compare the cost (counting only products and looking only at the leading term) of converting a skew-symmetric and a symmetric matrix to an upper-Hessenberg form using Householder reflections.

  comment

• 1.II.38C

  Numerical Analysis

  Part II, 2006

  (a) Define the Jacobi method with relaxation for solving the linear system $A x=b$.

  (b) Let $A$ be a symmetric positive definite matrix with diagonal part $D$ such that the matrix $2 D-A$ is also positive definite. Prove that the iteration always converges if the relaxation parameter $\omega$ is equal to 1 .

  (c) Let $A$ be the tridiagonal matrix with diagonal elements $a_{i i}=1$ and off-diagonal elements $a_{i+1, i}=a_{i, i+1}=1 / 4$. Prove that convergence occurs if $\omega$ satisfies $0<\omega \leqslant 4 / 3$. Explain briefly why the choice $\omega=1$ is optimal.

  [You may quote without proof any relevant result about the convergence of iterative methods and about the eigenvalues of matrices.]

  comment

• 2.II.38C

  Numerical Analysis

  Part II, 2006

  In the unit square the Poisson equation $\nabla^{2} u=f$, with zero Dirichlet boundary conditions, is being solved by the five-point formula using a square grid of mesh size $h=1 /(M+1)$,

  $$u_{i, j-1}+u_{i, j+1}+u_{i-1, j}+u_{i+1, j}-4 u_{i, j}=h^{2} f_{i, j} .$$

  Let $u(x, y)$ be the exact solution, and let $e_{i, j}=u_{i, j}-u(i h, j h)$ be the error of the five-point formula at the $(i, j)$ th grid point. Justifying each step, prove that

  $$\left[\sum_{i, j=1}^{M}\left|e_{i, j}\right|^{2}\right]^{1 / 2} \leqslant c h, \quad h \rightarrow 0$$

  where $c$ is some constant.

  comment

• 3.II.38C

  Numerical Analysis

  Part II, 2006

  (a) For the equation $y^{\prime}=f(t, y)$, consider the following multistep method with $s$ steps,

  $$\sum_{i=0}^{s} \rho_{i} y_{n+i}=h \sum_{i=0}^{s} \sigma_{i} f\left(t_{n+i}, y_{n+i}\right)$$

  where $h$ is the step size and $\rho_{i}, \sigma_{i}$ are specified constants with $\rho_{s}=1$. Prove that this method is of order $p$ if and only if

  $$\sum_{i=0}^{s} \rho_{i} Q\left(t_{n+i}\right)=h \sum_{i=0}^{s} \sigma_{i} Q^{\prime}\left(t_{n+i}\right)$$

  for any polynomial $Q$ of degree $\leqslant p$. Deduce that there is no $s$-step method of order $2 s+1$.

  [You may use the fact that, for any $a_{i}, b_{i}$, the Hermite interpolation problem

  $$Q\left(x_{i}\right)=a_{i}, \quad Q^{\prime}\left(x_{i}\right)=b_{i}, \quad i=0, \ldots, s$$

  is uniquely solvable in the space of polynomials of degree $2 s+1 .]$

  (b) State the Dahlquist equivalence theorem regarding the convergence of a multistep method. Determine all the values of the real parameter $a \neq 0$ for which the multistep method

  $$y_{n+3}+(2 a-3)\left[y_{n+2}-y_{n+1}\right]-y_{n}=h a\left[f_{n+2}+f_{n+1}\right]$$

  is convergent, and determine the order of convergence.

  comment

• 4.II.39C

  Numerical Analysis

  Part II, 2006

  The difference equation

  $$u_{m}^{n+1}=u_{m}^{n}+\frac{3}{2} \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)-\frac{1}{2} \mu\left(u_{m-1}^{n-1}-2 u_{m}^{n-1}+u_{m+1}^{n-1}\right),$$

  where $\mu=\Delta t /(\Delta x)^{2}$, is used to approximate a solution of the diffusion equation $u_{t}=u_{x x}$.

  (a) Prove that, as $\Delta t \rightarrow 0$ with constant $\mu$, the local error of the method is $\mathcal{O}(\Delta t)^{2}$.

  (b) Applying the Fourier stability test, show that the method is stable if and only if $\mu \leqslant \frac{1}{4}$.

  comment

• 1.II.38C

  Numerical Analysis

  Part II, 2007

  (a) For a numerical method to solve $y^{\prime}=f(t, y)$, define the linear stability domain and state when such a method is A-stable.

  (b) Determine all values of the real parameter $a$ for which the Runge-Kutta method

  $$\begin{aligned} k_{1} &=f\left(t_{n}+\left(\frac{1}{2}-a\right) h, y_{n}+h\left[\frac{1}{4} k_{1}+\left(\frac{1}{4}-a\right) k_{2}\right]\right) \\ k_{2} &=f\left(t_{n}+\left(\frac{1}{2}+a\right) h, y_{n}+h\left[\left(\frac{1}{4}+a\right) k_{1}+\frac{1}{4} k_{2}\right]\right) \\ y_{n+1} &=y_{n}+\frac{1}{2} h\left(k_{1}+k_{2}\right) \end{aligned}$$

  is A-stable.

  comment

• 2.II.38C

  Numerical Analysis

  Part II, 2007

  (a) State the Householder-John theorem and explain how it can be used to design iterative methods for solving a system of linear equations $A x=b$.

  (b) Let $A=L+D+U$ where $D$ is the diagonal part of $A$, and $L$ and $U$ are, respectively, the strictly lower and strictly upper triangular parts of $A$. Given a vector $b$, consider the following iterative scheme:

  $$(D+\omega L) x^{(k+1)}=(1-\omega) D x^{(k)}-\omega U x^{(k)}+\omega b$$

  Prove that if $A$ is a symmetric positive definite matrix, and $\omega \in(0,2)$, then the above iteration converges to the solution of the system $A x=b$.

  comment

• 3.II.38C

  Numerical Analysis

  Part II, 2007

  (a) Prove that all Toeplitz symmetric tridiagonal $M \times M$ matrices

  $$A=\left[\begin{array}{ccccc} a & b & 0 & \cdots & 0 \\ b & a & b & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & b & a & b \\ 0 & \cdots & 0 & b & a \end{array}\right]$$

  share the same eigenvectors $\left(v^{(k)}\right)_{k=1}^{M}$ with components $v_{i}^{(k)}=\sin \frac{k i \pi}{M+1}$, $i=1, \ldots, M$, and eigenvalues to be determined.

  (b) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0 \leqslant x \leqslant 1, \quad 0 \leqslant t \leqslant T$$

  is approximated by the Crank-Nicolson scheme

  $$\begin{aligned} u_{m}^{n+1}-\frac{1}{2} \mu\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right) &=u_{m}^{n}+\frac{1}{2} \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right), \\ \text { for } \quad m &=1, \ldots, M, \end{aligned}$$

  where $\mu=\Delta t /(\Delta x)^{2}, \quad \Delta x=1 /(M+1)$, and $u_{m}^{n}$ is an approximation to $u(m \Delta x, n \Delta t)$. Assuming that $u(0, t)=u(1, t)=0, \forall t, \quad$ show that the above scheme can be written in the form

  $$B u^{n+1}=C u^{n}, \quad 0 \leqslant n \leqslant(T / \Delta t)-1$$

  where $u^{n}=\left[u_{1}^{n}, \ldots, u_{M}^{n}\right]^{\top}$ and the real matrices $B$ and $C$ should be found. Using matrix analysis, find the range of $\mu$ for which the scheme is stable. [Do not use Fourier analysis.]

  comment

• 4.II.39C

  Numerical Analysis

  Part II, 2007

  (a) Suppose that $A$ is a real $n \times n$ matrix, and that $w \in \mathbb{R}^{n}$ and $\lambda_{1} \in \mathbb{R}$ are given so that $A w=\lambda_{1} w$. Further, let $S$ be a non-singular matrix such that $S w=c e^{(1)}$, where $e^{(1)}$ is the first coordinate vector and $c \neq 0$. Let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are $\lambda_{1}$ together with the eigenvalues of the bottom right $(n-1) \times(n-1)$ submatrix of $\widehat{A} .$

  (b) Suppose again that $A$ is a real $n \times n$ matrix, and that two linearly independent vectors $v, w \in \mathbb{R}^{n}$ are given such that the linear subspace $L\{v, w\}$ spanned by $v$ and $w$ is invariant under the action of $A$, i.e.,

  $$x \in L\{v, w\} \quad \Rightarrow \quad A x \in L\{v, w\}$$

  Denote by $V$ an $n \times 2$ matrix whose two columns are the vectors $v$ and $w$, and let $S$ be a non-singular matrix such that $R=S V$ is upper triangular, that is,

  $$R=S V=S \times\left[\begin{array}{cc} v_{1} & w_{1} \\ v_{2} & w_{2} \\ v_{3} & w_{3} \\ : & : \\ v_{n} & w_{n} \end{array}\right]=\left[\begin{array}{cc} r_{11} & r_{12} \\ 0 & r_{22} \\ 0 & 0 \\ : & : \\ 0 & 0 \end{array}\right]$$

  Again let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are the eigenvalues of the top left $2 \times 2$ submatrix of $\widehat{A}$together with the eigenvalues of the bottom right $(n-2) \times(n-2)$ submatrix of $\widehat{A}$

  comment

• 1.II.38C

  Numerical Analysis

  Part II, 2008

  The Poisson equation $\nabla^{2} u=f$ in the unit square $\Omega=[0,1] \times[0,1]$, with zero boundary conditions on $\partial \Omega$, is discretized with the nine-point formula

  $$\begin{aligned} \frac{10}{3} u_{m, n} &-\frac{2}{3}\left(u_{m+1, n}+u_{m-1, n}+u_{m, n+1}+u_{m, n-1}\right) \\ &-\frac{1}{6}\left(u_{m+1, n+1}+u_{m+1, n-1}+u_{m-1, n+1}+u_{m-1, n-1}\right)=-h^{2} f_{m, n}, \end{aligned}$$

  where $1 \leqslant m, n \leqslant M, u_{m, n} \approx u(m h, n h)$, and $(m h, n h)$ are grid points.

  (a) Prove that, for any ordering of the grid points, the method can be written as $A \mathbf{u}=\mathbf{b}$ with a symmetric positive-definite matrix $A$.

  (b) Describe the Jacobi method for solving a linear system of equations, and prove that it converges for the above system.

  [You may quote without proof the corollary of the Householder-John theorem regarding convergence of the Jacobi method.]

  comment

• 2.II.38C

  Numerical Analysis

  Part II, 2008

  The advection equation

  $$u_{t}=u_{x}, \quad x \in \mathbb{R}, \quad t \geqslant 0,$$

  is solved by the leapfrog scheme

  $$u_{m}^{n+1}=\mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right)+u_{m}^{n-1},$$

  where $n \geqslant 1$ and $\mu=\Delta t / \Delta x$ is the Courant number.

  (a) Determine the local error of the method.

  (b) Applying the Fourier technique, find the range of $\mu>0$ for which the method is stable.

  comment

• 3.II.38C

  Numerical Analysis

  Part II, 2008

  (a) A numerical method for solving the ordinary differential equation

  $$y^{\prime}(t)=f(t, y), \quad t \in[0, T], \quad y(0)=y_{0},$$

  generates for every $h>0$ a sequence $\left\{y_{n}\right\}$, where $y_{n}$ is an approximation to $y\left(t_{n}\right)$ and $t_{n}=n h$. Explain what is meant by the convergence of the method.

  (b) Prove from first principles that if the function $f$ is sufficiently smooth and satisfies the Lipschitz condition

  $$|f(t, x)-f(t, y)| \leqslant \lambda|x-y|, \quad x, y \in \mathbb{R}, \quad t \in[0, T],$$

  for some $\lambda>0$, then the trapezoidal rule

  $$y_{n+1}=y_{n}+\frac{1}{2} h\left[f\left(t_{n}, y_{n}\right)+f\left(t_{n+1}, y_{n+1}\right)\right]$$

  converges.

  comment

• 4.II.39C

  Numerical Analysis

  Part II, 2008

  Let $A \in \mathbb{R}^{n \times n}$ be a real matrix with $n$ linearly independent eigenvectors. When calculating eigenvalues of $A$, the sequence $\mathbf{x}^{(k)}, k=0,1,2, \ldots$, is generated by power method $\mathbf{x}^{(k+1)}=A \mathbf{x}^{(k)} /\left\|A \mathbf{x}^{(k)}\right\|$, where $\mathbf{x}^{(0)}$ is a real nonzero vector.

  (a) Describe the asymptotic properties of the sequence $\mathbf{x}^{(k)}$, both in the case where the eigenvalues $\lambda_{i}$ of $A$ satisfy $\left|\lambda_{i}\right|<\left|\lambda_{n}\right|, i=1, \ldots, n-1$, and in the case where $\left|\lambda_{i}\right|<\left|\lambda_{n-1}\right|=\left|\lambda_{n}\right|, i=1, \ldots, n-2$. In the latter case explain how the (possibly complexvalued) eigenvalues $\lambda_{n-1}, \lambda_{n}$ and their corresponding eigenvectors can be determined.

  (b) Let $n=3$, and suppose that, for a large $k$, we obtain the vectors

  $$\mathbf{y}_{k}=\mathbf{x}_{k}=\left[\begin{array}{l} 1 \\ 1 \\ 1 \end{array}\right], \quad \mathbf{y}_{k+1}=A \mathbf{x}_{k}=\left[\begin{array}{l} 2 \\ 3 \\ 4 \end{array}\right], \quad \mathbf{y}_{k+2}=A^{2} \mathbf{x}_{k}=\left[\begin{array}{l} 2 \\ 4 \\ 6 \end{array}\right]$$

  Find two eigenvalues of $A$ and their corresponding eigenvectors.

  comment

• Paper 3, Section II, B

  Numerical Analysis

  Part II, 2009

  Prove that all Toeplitz tridiagonal $M \times M$ matrices $A$ of the form

  $$A=\left[\begin{array}{rrrrr} a & b & & & \\ -b & a & b & & \\ & \ddots & \ddots & \ddots & \\ & & -b & a & b \\ & & & -b & a \end{array}\right]$$

  share the same eigenvectors $\left(\boldsymbol{v}^{(k)}\right)_{k=1}^{M}$, with the components $\boldsymbol{v}_{m}^{(k)}=i^{m} \sin \frac{k m \pi}{M+1}, m=$ $1, \ldots, M$, where $i=\sqrt{-1}$, and find their eigenvalues.

  The advection equation

  $$\frac{\partial u}{\partial t}=\frac{\partial u}{\partial x}, \quad 0 \leqslant x \leqslant 1, \quad 0 \leqslant t \leqslant T,$$

  is approximated by the Crank-Nicolson scheme

  $$u_{m}^{n+1}-u_{m}^{n}=\frac{1}{4} \mu\left(u_{m+1}^{n+1}-u_{m-1}^{n+1}\right)+\frac{1}{4} \mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right)$$

  where $\mu=\frac{\Delta t}{(\Delta x)^{2}}, \Delta x=\frac{1}{M+1}$, and $u_{m}^{n}$ is an approximation to $u(m \Delta x, n \Delta t)$. Assuming that $u(0, t)=u(1, t)=0$, show that the above scheme can be written in the form

  $$B \boldsymbol{u}^{n+1}=C \boldsymbol{u}^{n}, \quad 0 \leqslant n \leqslant T / \Delta t-1$$

  where $\boldsymbol{u}^{n}=\left[u_{1}^{n}, \ldots, u_{M}^{n}\right]^{T}$ and the real matrices $B$ and $C$ should be found. Using matrix analysis, find the range of $\mu$ for which the scheme is stable. [Fourier analysis is not acceptable.]

  comment

• Paper 2, Section II, B

  Numerical Analysis

  Part II, 2009

  The Poisson equation $\nabla^{2} u=f$ in the unit square $\Omega=[0,1] \times[0,1]$, equipped with appropriate boundary conditions on $\partial \Omega$, is discretized with the nine-point formula:

  $$\begin{aligned} \Gamma_{9}\left[u_{m, n}\right] &:=-\frac{10}{3} u_{m, n}+\frac{2}{3}\left(u_{m+1, n}+u_{m-1, n}+u_{m, n+1}+u_{m, n-1}\right) \\ &+\frac{1}{6}\left(u_{m+1, n+1}+u_{m+1, n-1}+u_{m-1, n+1}+u_{m-1, n-1}\right)=h^{2} f_{m, n} \end{aligned}$$

  where $1 \leqslant m, n \leqslant M, u_{m, n} \approx u(m h, n h)$, and $(m h, n h)$ are grid points.

  (i) Find the local error of approximation.

  (ii) Prove that the error is smaller if $f$ happens to satisfy the Laplace equation $\nabla^{2} f=0 .$

  (iii) Hence show that the modified nine-point scheme

  $$\begin{aligned} \Gamma_{9}\left[u_{m, n}\right] &=h^{2} f_{m, n}+\frac{1}{12} h^{2} \Gamma_{5}\left[f_{m, n}\right] \\ &:=h^{2} f_{m, n}+\frac{1}{12} h^{2}\left(f_{m+1, n}+f_{m-1, n}+f_{m, n+1}+f_{m, n-1}-4 f_{m, n}\right) \end{aligned}$$

  has the same smaller error as in (ii).

  [Hint. The nine-point discretization of $\nabla^{2} u$ can be written as

  $$\Gamma_{9}[u]=\left(\Gamma_{5}+\frac{1}{6} \Delta_{x}^{2} \Delta_{y}^{2}\right) u=\left(\Delta_{x}^{2}+\Delta_{y}^{2}+\frac{1}{6} \Delta_{x}^{2} \Delta_{y}^{2}\right) u$$

  where $\Gamma_{5}[u]=\left(\Delta_{x}^{2}+\Delta_{y}^{2}\right) u$ is the five-point discretization and

  $$\begin{aligned} &\Delta_{x}^{2} u(x, y):=u(x-h, y)-2 u(x, y)+u(x+h, y) \\ &\Delta_{y}^{2} u(x, y):=u(x, y-h)-2 u(x, y)+u(x, y+h) \end{aligned}$$

  comment

• Paper 1, Section II, B

  Numerical Analysis

  Part II, 2009

  (i) Define the Jacobi method with relaxation for solving the linear system $A \boldsymbol{x}=\boldsymbol{b}$.

  (ii) For $\boldsymbol{x}^{*}$ and $\boldsymbol{x}^{(\nu)}$ being the exact and the iterated solution, respectively, let $e^{(\nu)}:=\boldsymbol{x}^{(\nu)}-\boldsymbol{x}^{*}$ be the error and $H_{\omega}$ the iteration matrix, so that

  $$e^{(\nu+1)}=H_{\omega} e^{(\nu)}$$

  Express $H_{\omega}$ in terms of the matrix $A$, its diagonal part $D$ and the relaxation parameter $\omega$, and find the eigenvectors $\boldsymbol{v}_{k}$ and the eigenvalues $\lambda_{k}(\omega)$ of $H_{\omega}$ for the $n \times n$ tridiagonal matrix

  $$A=\left[\begin{array}{rrrrr} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{array}\right]$$

  [Hint: The eigenvectors and eigenvalues of $A$ are

  $$\left.\left(\boldsymbol{u}_{k}\right)_{i}=\sin \frac{k i \pi}{n+1}, \quad i=1, \ldots, n, \quad \lambda_{k}(A)=4 \sin ^{2} \frac{k \pi}{2(n+1)}, \quad k=1, \ldots, n .\right]$$

  (iii) For $A$ as above, let

  $$\boldsymbol{e}^{(\nu)}=\sum_{k=1}^{n} a_{k}^{(\nu)} \boldsymbol{v}_{k}$$

  be the expansion of the error with respect to the eigenvectors $\left(\boldsymbol{v}_{k}\right)$ of $H_{\omega}$.

  Find the range of parameter $\omega$ which provides convergence of the method for any $n$, and prove that, for any such $\omega$, the rate of convergence $e^{(\nu)} \rightarrow 0$ is not faster than $\left(1-c / n^{2}\right)^{\nu}$.

  (iv) Show that, for some $\omega$, the high frequency components $\left(\frac{n+1}{2} \leqslant k \leqslant n\right)$ of the error $e^{(\nu)}$ tend to zero much faster. Determine the optimal parameter $\omega_{*}$ which provides the largest suppression of the high frequency components per iteration, and find the corresponding attenuation factor $\mu_{*}$ (i.e. the least $\mu_{\omega}$ such that $\left|a_{k}^{(\nu+1)}\right| \leqslant \mu_{\omega}\left|a_{k}^{(\nu)}\right|$ for $\left.\frac{n+1}{2} \leqslant k \leqslant n\right)$.

  comment

• Paper 4, Section II, B

  Numerical Analysis

  Part II, 2009

  (a) For the $s$-step $s$-order Backward Differentiation Formula (BDF) for ordinary differential equations,

  $$\sum_{m=0}^{s} a_{m} y_{n+m}=h f_{n+s}$$

  express the polynomial $\rho(w)=\sum_{m=0}^{s} a_{m} w^{m}$ in a convenient explicit form.

  (b) Prove that the interval $(-\infty, 0)$ belongs to the linear stability domain of the 2-step BDF method.

  comment

• Paper 1, Section II, A

  Numerical Analysis

  Part II, 2010

  (a) State the Householder-John theorem and explain its relation to the convergence analysis of splitting methods for solving a system of linear equations $A x=b$ with a positive definite matrix $A$.

  (b) Describe the Jacobi method for solving a system $A x=b$, and deduce from the above theorem that if $A$ is a symmetric positive definite tridiagonal matrix,

  $$A=\left[\begin{array}{ccccc} a_{1} & c_{1} & & & \\ c_{1} & a_{2} & c_{2} & & \mathbf{0} \\ & \ddots & \ddots & \ddots & \\ \mathbf{0} & & c_{n-2} & a_{n-1} & c_{n-1} \\ & & & c_{n-1} & a_{n} \end{array}\right]$$

  then the Jacobi method converges.

  [Hint: At the last step, you may find it useful to consider two vectors $x=\left(x_{1}, x_{2}, \ldots, x_{n}\right)$ and $y=\left((-1) x_{1},(-1)^{2} x_{2}, \ldots,(-1)^{n} x_{n}\right)$.]

  comment

• Paper 2, Section II, A

  Numerical Analysis

  Part II, 2010

  The inverse discrete Fourier transform $\mathcal{F}_{n}^{-1}: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ is given by the formula

  $$\mathbf{x}=\mathcal{F}_{n}^{-1} \mathbf{y}, \quad \text { where } \quad x_{\ell}=\sum_{j=0}^{n-1} \omega_{n}^{j \ell} y_{j}, \quad \ell=0, \ldots, n-1$$

  Here, $\omega_{n}=\exp (2 \pi i / n)$ is the primitive root of unity of degree $n$, and $n=2^{p}, p=1,2, \ldots$

  (1) Show how to assemble $\mathbf{x}=\mathcal{F}_{2 m}^{-1} \mathbf{y}$ in a small number of operations if we already know the Fourier transforms of the even and odd portions of $\mathbf{y}$ :

  $$\mathbf{x}^{(\mathrm{E})}=\mathcal{F}_{m}^{-1} \mathbf{y}^{(\mathrm{E})}, \quad \mathbf{x}^{(\mathrm{O})}=\mathcal{F}_{m}^{-1} \mathbf{y}^{(\mathrm{O})}$$

  (2) Describe the Fast Fourier Transform (FFT) method for evaluating $\mathbf{x}$ and draw a relevant diagram for $n=8$.

  (3) Find the costs of the FFT for $n=2^{p}$ (only multiplications count).

  (4) For $n=4$, using the FFT technique, find $\mathbf{x}=\mathcal{F}_{4}^{-1} \mathbf{y}, \quad$ for $\quad \mathbf{y}=[1,1,-1,-1], \quad$ and $\quad \mathbf{y}=[1,-1,1,-1]$

  comment

• Paper 3, Section II, A

  Numerical Analysis

  Part II, 2010

  The Poisson equation $\nabla^{2} u=f$ in the unit square $\Omega=[0,1] \times[0,1], u=0$ on $\partial \Omega$, is discretized with the five-point formula

  $$u_{i, j-1}+u_{i, j+1}+u_{i+1, j}+u_{i-1, j}-4 u_{i, j}=h^{2} f_{i, j}$$

  where $1 \leqslant i, j \leqslant M, u_{i, j} \approx u(i h, j h)$ and $(i h, j h)$ are grid points.

  Let $u(x, y)$ be the exact solution, and let $e_{i, j}=u_{i, j}-u(i h, j h)$ be the error of the five-point formula at the $(i, j)$ th grid point. Justifying each step, prove that

  $$\|\mathbf{e}\|=\left[\sum_{i, j=1}^{M}\left|e_{i, j}\right|^{2}\right]^{1 / 2} \leqslant c h \quad \text { for sufficiently small } h>0,$$

  where $c$ is some constant independent of $h$.

  comment

• Paper 4, Section II, A

  Numerical Analysis

  Part II, 2010

  An $s$-stage explicit Runge-Kutta method of order $p$, with constant step size $h>0$, is applied to the differential equation $y^{\prime}=\lambda y, t \geqslant 0$.

  (a) Prove that

  $$y_{n+1}=P_{s}(\lambda h) y_{n}$$

  where $P_{s}$ is a polynomial of degree $s$.

  (b) Prove that the order $p$ of any $s$-stage explicit Runge-Kutta method satisfies the inequality $p \leqslant s$ and, for $p=s$, write down an explicit expression for $P_{s}$.

  (c) Prove that no explicit Runge-Kutta method can be A-stable.

  comment

• Paper 1, Section II, A

  Numerical Analysis

  Part II, 2011

  The nine-point method for the Poisson equation $\nabla^{2} u=f$ (with zero Dirichlet boundary conditions) in a square, reads

  $$\begin{array}{r} \frac{2}{3}\left(u_{i-1, j}+u_{i+1, j}+u_{i, j-1}+u_{i, j+1}\right)+\frac{1}{6}\left(u_{i-1, j-1}+u_{i-1, j+1}+u_{i+1, j-1}+u_{i+1, j+1}\right) \\ -\frac{10}{3} u_{i, j}=h^{2} f_{i, j}, \quad i, j=1, \ldots, m, \end{array}$$

  where $u_{0, j}=u_{m+1, j}=u_{i, 0}=u_{i, m+1}=0$, for all $i, j=0, \ldots, m+1$.

  (i) By arranging the two-dimensional arrays $\left\{u_{i, j}\right\}_{i, j=1, \ldots, m}$ and $\left\{f_{i, j}\right\}_{i, j=1, \ldots, m}$ into column vectors $u \in \mathbb{R}^{m^{2}}$ and $b \in \mathbb{R}^{m^{2}}$ respectively, the linear system above takes the matrix form $A u=b$. Prove that, regardless of the ordering of the points on the grid, the matrix $A$ is symmetric and negative definite.

  (ii) Formulate the Jacobi method with relaxation for solving the above linear system.

  (iii) Prove that the iteration converges if the relaxation parameter $\omega$ is equal to $1 .$

  [You may quote without proof any relevant result about convergence of iterative methods.]

  comment

• Paper 2, Section II, A

  Numerical Analysis

  Part II, 2011

  Let $A \in \mathbb{R}^{n \times n}$ be a real matrix with $n$ linearly independent eigenvectors. The eigenvalues of $A$ can be calculated from the sequence $x^{(k)}, k=0,1, \ldots$, which is generated by the power method

  $$x^{(k+1)}=\frac{A x^{(k)}}{\left\|A x^{(k)}\right\|},$$

  where $x^{(0)}$ is a real nonzero vector.

  (i) Describe the asymptotic properties of the sequence $x^{(k)}$ in the case that the eigenvalues $\lambda_{i}$ of $A$ satisfy $\left|\lambda_{i}\right|<\left|\lambda_{n}\right|, i=1, \ldots, n-1$, and the eigenvectors are of unit length.

  (ii) Present the implementation details for the power method for the setting in (i) and define the Rayleigh quotient.

  (iii) Let $A$ be the $3 \times 3$ matrix

  $$A=\lambda I+P, \quad P=\left(\begin{array}{lll} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right)$$

  where $\lambda$ is real and nonzero. Find an explicit expression for $A^{k}, k=1,2,3, \ldots$

  Let the sequence $x^{(k)}$ be generated by the power method as above. Deduce from your expression for $A^{k}$ that the first and second components of $x^{(k+1)}$ tend to zero as $k \rightarrow \infty$. Further show that this implies $A x^{(k+1)}-\lambda x^{(k+1)} \rightarrow 0$ as $k \rightarrow \infty$.

  comment

• Paper 3, Section II, A

  Numerical Analysis

  Part II, 2011

  (i) The difference equation

  $$u_{i}^{n+1}=u_{i}^{n}+\frac{3}{2} \mu\left(u_{i-1}^{n}-2 u_{i}^{n}+u_{i+1}^{n}\right)-\frac{1}{2} \mu\left(u_{i-1}^{n-1}-2 u_{i}^{n-1}+u_{i+1}^{n-1}\right)$$

  where $\mu=\Delta t /(\Delta x)^{2}$, is the basic equation used in the second-order AdamsBashforth method and can be employed to approximate a solution of the diffusion equation $u_{t}=u_{x x}$. Prove that, as $\Delta t \rightarrow 0$ with constant $\mu$, the local error of the method is $O(\Delta t)^{2}$.

  (ii) By applying the Fourier stability test, show that the above method is stable if and only if $\mu \leqslant 1 / 4$.

  (iii) Define the leapfrog scheme to approximate the diffusion equation and prove that it is unstable for every choice of $\mu>0$.

  comment

• Paper 4, Section II, A

  Numerical Analysis

  Part II, 2011

  (i) Consider the Poisson equation

  $$\nabla^{2} u=f, \quad-1 \leqslant x, y \leqslant 1$$

  with the periodic boundary conditions

  and the normalization condition

  $$\int_{-1}^{1} \int_{-1}^{1} u(x, y) d x d y=0$$

  Moreover, $f$ is analytic and obeys the periodic boundary conditions $f(-1, y)=$ $f(1, y), f(x,-1)=f(x, 1),-1 \leqslant x, y \leqslant 1 .$

  Derive an explicit expression of the approximation of a solution $u$ by means of a spectral method. Explain the term convergence with spectral speed and state its validity for the approximation of $u$.

  (ii) Consider the second-order linear elliptic partial differential equation

  $$\nabla \cdot(a \nabla u)=f, \quad-1 \leqslant x, y \leqslant 1$$

  with the periodic boundary conditions and normalization condition specified in (i). Moreover, $a$ and $f$ are given by

  $$a(x, y)=\cos (\pi x)+\cos (\pi y)+3, \quad f(x, y)=\sin (\pi x)+\sin (\pi y)$$

  [Note that $a$ is a positive analytic periodic function.]

  Construct explicitly the linear algebraic system that arises from the implementation of a spectral method to the above equation.

  $$\begin{aligned} & u(-1, y)=u(1, y), \quad u_{x}(-1, y)=u_{x}(1, y), \quad-1 \leqslant y \leqslant 1, \\ & u(x,-1)=u(x, 1), \quad u_{y}(x,-1)=u_{y}(x, 1), \quad-1 \leqslant x \leqslant 1 \end{aligned}$$

  comment

• Paper 4, Section II, $39 \mathrm{D}$

  Numerical Analysis

  Part II, 2012

  (i) Formulate the conjugate gradient method for the solution of a system $A \mathbf{x}=\mathbf{b}$ with $A \in \mathbb{R}^{n \times n}$ and $\mathbf{b} \in \mathbb{R}^{n}, n>0$.

  (ii) Prove that if the conjugate gradient method is applied in exact arithmetic then, for any $\mathbf{x}^{(0)} \in \mathbb{R}^{n}$, termination occurs after at most $n$ iterations.

  (iii) The polynomial $p(x)=x^{m}+\sum_{i=0}^{m-1} c_{i} x^{i}$ is the minimal polynomial of the $n \times n$ matrix $A$ if it is the polynomial of lowest degree that satisfies $p(A)=0 . \quad[$ Note that $m \leqslant n$.] Give an example of a $3 \times 3$ symmetric positive definite matrix with a quadratic minimal polynomial.

  Prove that (in exact arithmetic) the conjugate gradient method requires at most $m$ iterations to calculate the exact solution of $A \mathbf{x}=\mathbf{b}$, where $m$ is the degree of the minimal polynomial of $A$.

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part II, 2012

  (i) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0 \leqslant x \leqslant 1, t \geqslant 0$$

  with the initial condition $u(x, 0)=\phi(x), 0 \leqslant x \leqslant 1$, and with zero boundary conditions at $x=0$ and $x=1$, can be solved numerically by the method

  $$u_{m}^{n+1}=u_{m}^{n}+\mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right), \quad m=1,2, \ldots, M, n \geqslant 0$$

  where $\Delta x=1 /(M+1), \mu=\Delta t /(\Delta x)^{2}$, and $u_{m}^{n} \approx u(m \Delta x, n \Delta t)$. Prove that $\mu \leqslant 1 / 2$ implies convergence.

  (ii) By discretising the diffusion equation and employing the same notation as in (i) above, determine [without using Fourier analysis] conditions on $\mu$ and the constant $\alpha$ such that the method

  $$u_{m}^{n+1}-\frac{1}{2}(\mu-\alpha)\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right)=u_{m}^{n}+\frac{1}{2}(\mu+\alpha)\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)$$

  is stable.

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part II, 2012

  The inverse discrete Fourier transform $\mathcal{F}_{n}^{-1}: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ is given by the formula

  $$\mathbf{x}=\mathcal{F}_{n}^{-1} \mathbf{y}, \quad \text { where } \quad x_{l}=\sum_{j=0}^{n-1} \omega_{n}^{j l} y_{j}, \quad l=0, \ldots, n-1$$

  Here, $\omega_{n}=\exp (2 \pi i / n)$ is the primitive root of unity of degree $n$ and $n=2^{p}, p=1,2, \ldots$

  (i) Show how to assemble $\mathbf{x}=\mathcal{F}_{2 m}^{-1} \mathbf{y}$ in a small number of operations if the Fourier transforms of the even and odd parts of $\mathbf{y}$,

  $$\mathbf{x}^{(E)}=\mathcal{F}_{m}^{-1} \mathbf{y}^{(E)}, \quad \mathbf{x}^{(O)}=\mathcal{F}_{m}^{-1} \mathbf{y}^{(O)}$$

  are already known.

  (ii) Describe the Fast Fourier Transform (FFT) method for evaluating $\mathbf{x}$, and draw a relevant diagram for $n=8$.

  (iii) Find the costs of the FFT method for $n=2^{p}$ (only multiplications count).

  (iv) For $n=4$ use the FFT method to find $\mathbf{x}=\mathcal{F}_{4}^{-1} \mathbf{y}$ when:

  (a) $\mathbf{y}=(1,1,-1,-1)$,

  (b) $\mathbf{y}=(1,-1,1,-1)$.

  comment

• Paper 1, Section II, D

  Numerical Analysis

  Part II, 2012

  The Poisson equation $u_{x x}=f$ in the unit interval $\Omega=[0,1], u=0$ on $\partial \Omega$ is discretised with the formula

  $$u_{i-1}+u_{i+1}-2 u_{i}=h^{2} f_{i}$$

  where $1 \leqslant i \leqslant n, u_{i} \approx u(i h)$ and $i h$ are the grid points.

  (i) Define the above system of equations in vector form $A \mathbf{u}=\mathbf{b}$ and describe the relaxed Jacobi method with relaxation parameter $\omega$ for solving this linear system. For $\mathbf{x}^{*}$ and $\mathbf{x}^{(\nu)}$ being the exact solution and the iterated solution respectively, let $\mathbf{e}^{(\nu)}=\mathbf{x}^{(\nu)}-\mathbf{x}^{*}$ be the error and $H_{\omega}$ the iteration matrix, so that

  $$\mathbf{e}^{(\nu+1)}=H_{\omega} \mathbf{e}^{(\nu)}$$

  Express $H_{\omega}$ in terms of the matrix $A$, the diagonal part $D$ of $A$ and $\omega$, and find the eigenvectors $\mathbf{v}_{k}$ and the eigenvalues $\lambda_{k}(\omega)$ of $H_{\omega}$.

  (ii) For $A$ as above, let

  $$\mathbf{e}^{(\nu)}=\sum_{k=1}^{n} a_{k}^{(\nu)} \mathbf{v}_{k}$$

  be the expansion of the error with respect to the eigenvectors of $H_{\omega}$. Derive conditions on $\omega$ such that the method converges for any $n$, and prove that, for any such $\omega$, the rate of convergence of $\mathbf{e}^{(\nu)} \rightarrow 0$ is not faster than $\left(1-c / n^{2}\right)^{\nu}$.

  (iii) Show that, for some $\omega$, the high frequency components $\left(\frac{n+1}{2} \leqslant k \leqslant n\right)$ of the error $\mathbf{e}^{(\nu)}$ tend to zero much faster than $\left(1-c / n^{2}\right)^{\nu}$. Determine the optimal parameter $\omega_{*}$ which provides the largest suppression of the high frequency components per iteration, and find the corresponding attenuation factor $\mu_{*}$ (i.e., the least $\mu_{\omega}$ such that $\left|a_{k}^{(\nu+1)}\right| \leqslant \mu_{\omega}\left|a_{k}^{(\nu)}\right|$ for $\left.\frac{n+1}{2} \leqslant k \leqslant n\right)$.

  comment

• Paper 4, Section II, C

  Numerical Analysis

  Part II, 2013

  Consider the solution of the two-point boundary value problem

  $$(2-\sin \pi x) u^{\prime \prime}+u=1, \quad-1 \leqslant x \leqslant 1$$

  with periodic boundary conditions at $x=-1$ and $x=1$. Construct explicitly the linear algebraic system that arises from the application of a spectral method to the above equation.

  The Fourier coefficients of $u$ are defined by

  $$\hat{u}_{n}=\frac{1}{2} \int_{-1}^{1} u(\tau) e^{-i \pi n \tau} d \tau$$

  Prove that the computation of the Fourier coefficients for the truncated system with $-N / 2+1 \leqslant n \leqslant N / 2$ (where $N$ is an even and positive integer, and assuming that $\hat{u}_{n}=0$ outside this range of $n$ ) reduces to the solution of a tridiagonal system of algebraic equations, which you should specify.

  Explain the term convergence with spectral speed and justify its validity for the derived approximation of $u$.

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part II, 2013

  Consider the advection equation $u_{t}=u_{x}$ on the unit interval $x \in[0,1]$ and $t \geqslant 0$, where $u=u(x, t)$, subject to the initial condition $u(x, 0)=\varphi(x)$ and the boundary condition $u(1, t)=0$, where $\varphi$ is a given smooth function on $[0,1]$.

  (i) We commence by discretising the advection equation above with finite differences on the equidistant space-time grid $\{(m \Delta x, n \Delta t), m=0, \ldots, M+1, n=0, \ldots, T\}$ with $\Delta x=1 /(M+1)$ and $\Delta t>0$. We obtain an equation for $u_{m}^{n} \approx u(m \Delta x, n \Delta t)$ that reads

  $$u_{m}^{n+1}=u_{m}^{n}+\frac{1}{2} \mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right), \quad m=1, \ldots, M, n \in \mathbb{Z}^{+}$$

  with the condition $u_{0}^{n}=0$ for all $n \in \mathbb{Z}^{+}$and $\mu=\Delta t / \Delta x$.

  What is the order of approximation (that is, the order of the local error) in space and time of the above discrete solution to the exact solution of the advection equation? Write the scheme in matrix form and deduce for which choices of $\mu$ this approximation converges to the exact solution. State (without proof) any theorems you use. [You may use the fact that for a tridiagonal $M \times M$ matrix

  $$\left(\begin{array}{cccc} \alpha & \beta & 0 & 0 \\ -\beta & \ddots & \ddots & 0 \\ 0 & \ddots & \ddots & \beta \\ 0 & 0 & -\beta & \alpha \end{array}\right)$$

  the eigenvalues are given by $\lambda_{\ell}=\alpha+2 i \beta \cos \frac{\ell \pi}{M+1}$.]

  (ii) How does the order change when we replace the central difference approximation of the first derivative in space by forward differences, that is $u_{m+1}^{n}-u_{m}^{n}$ instead of $\left(u_{m+1}^{n}-u_{m-1}^{n}\right) / 2 ?$ For which choices of $\mu$ is this new scheme convergent?

  (iii) Instead of the approximation in (i) we consider the following method for numerically solving the advection equation,

  $$u_{m}^{n+1}=\mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right)+u_{m}^{n-1}$$

  where we additionally assume that $u_{m}^{1}$ is given. What is the order of this method for a fixed $\mu$ ?

  comment

• Paper 3, Section II, C

  Numerical Analysis

  Part II, 2013

  (i) Suppose that $A$ is a real $n \times n$ matrix, and that $\mathbf{w} \in \mathbb{R}^{n}$ and $\lambda_{1} \in \mathbb{R}$ are given so that $A \mathbf{w}=\lambda_{1} \mathbf{w}$. Further, let $S$ be a non-singular matrix such that $S \mathbf{w}=c \mathbf{e}_{1}$, where $\mathbf{e}_{1}$ is the first coordinate vector and $c \neq 0$. Let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are $\lambda_{1}$ together with the eigenvalues of the bottom right $(n-1) \times(n-1)$ submatrix of $\widehat{A}$.

  (ii) Suppose again that $A$ is a real $n \times n$ matrix, and that two linearly independent vectors $\mathbf{v}, \mathbf{w} \in \mathbb{R}^{n}$ are given such that the linear subspace $\mathcal{L}\{\mathbf{v}, \mathbf{w}\}$ spanned by $\mathbf{v}$ and $w$ is invariant under the action of $A$, that is

  $$x \in \mathcal{L}\{\mathbf{v}, \mathbf{w}\} \quad \Rightarrow \quad A x \in \mathcal{L}\{\mathbf{v}, \mathbf{w}\} .$$

  Denote by $V$ an $n \times 2$ matrix whose two columns are the vectors $\mathbf{v}$ and $\mathbf{w}$, and let $S$ be a non-singular matrix such that $R=S V$ is upper triangular, that is

  $$R=S V=S \times\left(\begin{array}{cc} v_{1} & w_{1} \\ v_{2} & w_{2} \\ \vdots & \vdots \\ v_{n} & w_{n} \end{array}\right)=\left(\begin{array}{cc} r_{11} & r_{12} \\ 0 & r_{22} \\ 0 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{array}\right)$$

  Again, let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are the eigenvalues of the top left $2 \times 2$ submatrix of $\widehat{A}$together with the eigenvalues of the bottom right $(n-2) \times(n-2)$ submatrix of $\widehat{A} .$

  comment

• Paper 1, Section II, 40C

  Numerical Analysis

  Part II, 2013

  Let

  $$A(\alpha)=\left(\begin{array}{ccc} 1 & \alpha & \alpha \\ \alpha & 1 & \alpha \\ \alpha & \alpha & 1 \end{array}\right), \quad \alpha \in \mathbb{R}$$

  (i) For which values of $\alpha$ is $A(\alpha)$ positive definite?

  (ii) Formulate the Gauss-Seidel method for the solution $\mathbf{x} \in \mathbb{R}^{3}$ of a system

  $$A(\alpha) \mathbf{x}=\mathbf{b}$$

  with $A(\alpha)$ as defined above and $\mathbf{b} \in \mathbb{R}^{3}$. Prove that the Gauss-Seidel method converges to the solution of the above system whenever $A$ is positive definite. [You may state and use the Householder-John theorem without proof.]

  (iii) For which values of $\alpha$ does the Jacobi iteration applied to the solution of the above system converge?

  comment

• Paper 4, Section II, D

  Numerical Analysis

  Part II, 2014

  Let $A$ be a real symmetric $n \times n$ matrix with $n$ distinct real eigenvalues $\lambda_{1}<\lambda_{2}<$ $\cdots<\lambda_{n}$ and a corresponding orthogonal basis of normalized real eigenvectors $\left\{\mathbf{w}_{i}\right\}_{i=1}^{n}$.

  (i) Let $s \in \mathbb{R}$ satisfy $s<\lambda_{1}$. Given a unit vector $\mathbf{x}^{(0)} \in \mathbb{R}^{n}$, the iteration scheme

  $$\begin{gathered} (A-s I) \mathbf{y}=\mathbf{x}^{(k)} \\ \mathbf{x}^{(k+1)}=\mathbf{y} /\|\mathbf{y}\| \end{gathered}$$

  generates a sequence of vectors $\mathbf{x}^{(k+1)}$ for $k=0,1,2, \ldots$. Assuming that $\mathbf{x}^{(0)}=\sum c_{i} \mathbf{w}_{i}$ with $c_{1} \neq 0$, prove that $\mathbf{x}^{(k)}$ tends to $\pm \mathbf{w}_{1}$ as $k \rightarrow \infty$. What happens to $\mathbf{x}^{(k)}$ if $s>\lambda_{1}$ ? [Consider all cases.]

  (ii) Describe how to implement an inverse-iteration algorithm to compute the eigenvalues and eigenvectors of $A$, given some initial estimates for the eigenvalues.

  (iii) Let $n=2$. For iterates $\mathbf{x}^{(k)}$ of an inverse-iteration algorithm with a fixed value of $s \neq \lambda_{1}, \lambda_{2}$, show that if

  $$\mathbf{x}^{(k)}=\left(\mathbf{w}_{1}+\epsilon_{k} \mathbf{w}_{2}\right) /\left(1+\epsilon_{k}^{2}\right)^{1 / 2}$$

  where $\left|\epsilon_{k}\right|$ is small, then $\left|\epsilon_{k+1}\right|$ is of the same order of magnitude as $\left|\epsilon_{k}\right|$.

  (iv) Let $n=2$ still. Consider the iteration scheme

  $$s_{k}=\left(\mathbf{x}^{(k)}, A \mathbf{x}^{(k)}\right), \quad\left(A-s_{k} I\right) \mathbf{y}=\mathbf{x}^{(k)}, \quad \mathbf{x}^{(k+1)}=\mathbf{y} /\|\mathbf{y}\|$$

  for $k=0,1,2, \ldots$, where $(,$, denotes the inner product. Show that with this scheme $\left|\epsilon_{k+1}\right|=\left|\epsilon_{k}\right|^{3} .$

  comment

• Paper 2, Section II, D

  Numerical Analysis

  Part II, 2014

  Consider the one-dimensional advection equation

  $$u_{t}=u_{x}, \quad-\infty<x<\infty, \quad t \geqslant 0$$

  subject to an initial condition $u(x, 0)=\varphi(x)$. Consider discretization of this equation with finite differences on an equidistant space-time $\left\{(m h, n k), m \in \mathbb{Z}, n \in \mathbb{Z}^{+}\right\}$with step size $h>0$ in space and step size $k>0$ in time. Define the Courant number $\mu$ and explain briefly how such a discretization can be used to derive numerical schemes in which solutions $u_{m}^{n} \approx u(m h, n k), m \in \mathbb{Z}$ and $n \in \mathbb{Z}^{+}$satisfy equations of the form

  $$\sum_{i=r}^{s} a_{i} u_{m+i}^{n+1}=\sum_{i=r}^{s} b_{i} u_{m+i}^{n}$$

  where the coefficients $a_{i}, b_{i}$ are independent of $m, n$.

  (i) Define the order of a numerical scheme such as (1). Define what a convergent numerical scheme is. Explain the notion of stability and state the Lax equivalence theorem that connects convergence and stability of numerical schemes for linear partial differential equations.

  (ii) Consider the following example of (1):

  $$u_{m}^{n+1}=u_{m}^{n}+\frac{\mu}{2}\left(u_{m+1}^{n}-u_{m-1}^{n}\right)+\frac{\mu^{2}}{2}\left(u_{m+1}^{n}-2 u_{m}^{n}+u_{m-1}^{n}\right) .$$

  Determine conditions on $\mu$ such that the scheme $(2)$ is stable and convergent. What is the order of this scheme?

  comment

• Paper 3, Section II, D

  Numerical Analysis

  Part II, 2014

  Consider the linear system

  $$A x=b$$

  where $A \in \mathbb{R}^{n \times n}$ and $b, x \in \mathbb{R}^{n}$.

  (i) Define the Jacobi iteration method with relaxation parameter $\omega$ for solving (1).

  (ii) Assume that $A$ is a symmetric positive-definite matrix whose diagonal part $D$ is such that the matrix $2 D-A$ is also positive definite. Prove that the relaxed Jacobi iteration method always converges if the relaxation parameter $\omega$ is equal to $1 .$

  (iii) Let $A$ be the tridiagonal matrix with diagonal elements $a_{i i}=\alpha$ and off-diagonal elements $a_{i+1, i}=a_{i, i+1}=\beta$, where $0<\beta<\frac{1}{2} \alpha$. For which values of $\omega$ (expressed in terms of $\alpha$ and $\beta$ ) does the relaxed Jacobi iteration method converge? What choice of $\omega$ gives the optimal convergence speed?

  [You may quote without proof any relevant results about the convergence of iterative methods and about the eigenvalues of matrices.]

  comment

• Paper 1, Section II, 40D

  Numerical Analysis

  Part II, 2014

  (i) Consider the numerical approximation of the boundary-value problem

  $$\begin{gathered} u^{\prime \prime}=f, \quad u:[0,1] \rightarrow \mathbb{R}, \\ u(0)=\varphi_{0}, \quad u(1)=\varphi_{1}, \end{gathered}$$

  where $\varphi_{0}, \varphi_{1}$ are given constants and $f$ is a given smooth function on $[0,1]$. A grid $\left\{x_{1}, x_{2}, \ldots, x_{N}\right\}, N \geqslant 3$, on $[0,1]$ is given by

  $$x_{1}=\alpha_{1} h, \quad x_{i}=x_{i-1}+h \text { for } i=2, \ldots, N-1, \quad x_{N}=1-\alpha_{2} h$$

  where $0<\alpha_{1}, \alpha_{2}<1, \alpha_{1}+\alpha_{2}=1$ and $h=1 / N$. Derive finite-difference approximations for $u^{\prime \prime}\left(x_{i}\right)$, for $i=1, \ldots, N$, using at most one neighbouring grid point of $x_{i}$ on each side. Hence write down a numerical scheme to solve the problem, displaying explicitly the entries of the system matrix $A$ in the resulting system of linear equations $A u=b$, $A \in \mathbb{R}^{N \times N}, u, b \in \mathbb{R}^{N}$. What is the overall order of this numerical scheme? Explain briefly one strategy by which the order could be improved with the same grid.

  (ii) Consider the numerical approximation of the boundary-value problem

  $$\begin{aligned} &\nabla^{2} u=f, \quad u: \Omega \rightarrow \mathbb{R}, \\ &u(x)=0 \text { for all } x \in \partial \Omega \end{aligned}$$

  where $\Omega \subset \mathbb{R}^{2}$ is an arbitrary, simply connected bounded domain with smooth boundary $\partial \Omega$, and $f$ is a given smooth function. Define the 9-point formula used to approximate the Laplacian. Using this formula and an equidistant grid inside $\Omega$, define a numerical scheme for which the system matrix is symmetric and negative definite. Prove that the system matrix of your scheme has these properties for all choices of ordering of the grid points.

  Part II, $2014 \quad$ List of Questions

  [TURN OVER

  comment

• Paper 4, Section II, E

  Numerical Analysis

  Part II, 2015

  (a) Define the $m$ th Krylov space $K_{m}(A, v)$ for $A \in \mathbb{R}^{n \times n}$ and $0 \neq v \in \mathbb{R}^{n}$. Letting $\delta_{m}$ be the dimension of $K_{m}(A, v)$, prove the following results.

  (i) There exists a positive integer $s \leqslant n$ such that $\delta_{m}=m$ for $m \leqslant s$ and $\delta_{m}=s$ for $m>s$.

  (ii) If $v=\sum_{i=1}^{s^{\prime}} c_{i} w_{i}$, where $w_{i}$ are eigenvectors of $A$ for distinct eigenvalues and all $c_{i}$ are nonzero, then $s=s^{\prime}$.

  (b) Define the term residual in the conjugate gradient (CG) method for solving a system $A x=b$ with symmetric positive definite $A$. Explain (without proof) the connection to Krylov spaces and prove that for any right-hand side $b$ the CG method finds an exact solution after at most $t$ steps, where $t$ is the number of distinct eigenvalues of $A$. [You may use without proof known properties of the iterates of the CG method.]

  Define what is meant by preconditioning, and explain two ways in which preconditioning can speed up convergence. Can we choose the preconditioner so that the CG method requires only one step? If yes, is it a reasonable method for speeding up the computation?

  comment

• Paper 2, Section II, E

  Numerical Analysis

  Part II, 2015

  (a) The boundary value problem $-\Delta u+c u=f$ on the unit square $[0,1]^{2}$ with zero boundary conditions and scalar constant $c>0$ is discretised using finite differences as

  $$\begin{array}{r} -u_{i-1, j}-u_{i+1, j}-u_{i, j-1}-u_{i, j+1}+4 u_{i, j}+c h^{2} u_{i, j}=h^{2} f(i h, j h), \\ i, j=1, \ldots, m \end{array}$$

  with $h=1 /(m+1)$. Show that for the resulting system $A u=b$, for a suitable matrix $A$ and vectors $u$ and $b$, both the Jacobi and Gauss-Seidel methods converge. [You may cite and use known results on the discretised Laplace operator and on the convergence of iterative methods.]

  Define the Jacobi method with relaxation parameter $\omega$. Find the eigenvalues $\lambda_{k, l}$ of the iteration matrix $H_{\omega}$ for the above problem and show that, in order to ensure convergence for all $h$, the condition $0<\omega \leqslant 1$ is necessary.

  [Hint: The eigenvalues of the discretised Laplace operator in two dimensions are $4\left(\sin ^{2} \frac{\pi k h}{2}+\sin ^{2} \frac{\pi l h}{2}\right)$ for integers $k, l$.]

  (b) Explain the components and steps in a multigrid method for solving the Poisson equation, discretised as $A_{h} u_{h}=b_{h}$. If we use the relaxed Jacobi method within the multigrid method, is it necessary to choose $\omega \neq 1$ to get fast convergence? Explain why or why not.

  comment

• Paper 3, Section II, E

  Numerical Analysis

  Part II, 2015

  (a) Given the finite-difference recurrence

  $$\sum_{k=r}^{s} a_{k} u_{m+k}^{n+1}=\sum_{k=r}^{s} b_{k} u_{m+k}^{n}, \quad m \in \mathbb{Z}, n \in \mathbb{Z}^{+}$$

  that discretises a Cauchy problem, the amplification factor is defined by

  $$H(\theta)=\left(\sum_{k=r}^{s} b_{k} e^{i k \theta}\right) /\left(\sum_{k=r}^{s} a_{k} e^{i k \theta}\right)$$

  Show how $H(\theta)$ acts on the Fourier transform $\hat{u}^{n}$ of $u^{n}$. Hence prove that the method is stable if and only if $|H(\theta)| \leqslant 1$ for all $\theta \in[-\pi, \pi]$.

  (b) The two-dimensional diffusion equation

  $$u_{t}=u_{x x}+c u_{y y}$$

  for some scalar constant $c>0$ is discretised with the forward Euler scheme

  $$u_{i, j}^{n+1}=u_{i, j}^{n}+\mu\left(u_{i+1, j}^{n}-2 u_{i, j}^{n}+u_{i-1, j}^{n}+c u_{i, j+1}^{n}-2 c u_{i, j}^{n}+c u_{i, j-1}^{n}\right)$$

  Using Fourier stability analysis, find the range of values $\mu>0$ for which the scheme is stable.

  comment

• Paper 1, Section II, 38E

  Numerical Analysis

  Part II, 2015

  (a) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial}{\partial x}\left(a(x) \frac{\partial u}{\partial x}\right) \quad \text { in } \quad 0 \leqslant x \leqslant 1, \quad t \geqslant 0$$

  with the initial condition $u(x, 0)=\phi(x)$ in $0 \leqslant x \leqslant 1$ and zero boundary conditions at $x=0$ and $x=1$, is solved by the finite-difference method

  $$\begin{array}{r} u_{m}^{n+1}=u_{m}^{n}+\mu\left[a_{m-\frac{1}{2}} u_{m-1}^{n}-\left(a_{m-\frac{1}{2}}+a_{m+\frac{1}{2}}\right) u_{m}^{n}+a_{m+\frac{1}{2}} u_{m+1}^{n}\right] \\ m=1,2, \ldots, M \end{array}$$

  where $\mu=k / h^{2}, k=\Delta t, h=1 /(M+1), u_{m}^{n} \approx u(m h, n k)$, and $a_{\alpha}=a(\alpha h)$.

  Assuming that the function $a$ and the exact solution are sufficiently smooth, prove that the exact solution satisfies the numerical scheme with error $O\left(h^{3}\right)$ for constant $\mu$.

  (b) For the problem in part (a), assume that there exist $0<a_{-}<a_{+}<\infty$ such that $a_{-} \leqslant a(x) \leqslant a_{+}$for all $0 \leqslant x \leqslant 1$. State (without proof) the Gershgorin theorem and prove that the method is stable for $0<\mu \leqslant 1 /\left(2 a_{+}\right)$.

  comment

• Paper 4, Section II, B

  Numerical Analysis

  Part II, 2016

  (a) Describe an implementation of the power method for determining the eigenvalue of largest modulus and its associated eigenvector for a matrix that has a unique eigenvalue of largest modulus.

  Now let $A$ be a real $n \times n$ matrix with distinct eigenvalues satisfying $\left|\lambda_{n}\right|=\left|\lambda_{n-1}\right|$ and $\left|\lambda_{n}\right|>\left|\lambda_{i}\right|, i=1, \ldots, n-2$. The power method is applied to $A$, with an initial condition $\mathbf{x}^{(0)}=\sum_{i=1}^{n} c_{i} \mathbf{w}_{i}$ such that $c_{n-1} c_{n} \neq 0$, where $\mathbf{w}_{i}$ is the eigenvector associated with $\lambda_{i}$. Show that the power method does not converge. Explain why $\mathbf{x}^{(k)}, \mathbf{x}^{(k+1)}$ and $\mathbf{x}^{(k+2)}$ become linearly dependent as $k \rightarrow \infty$.

  (b) Consider the following variant of the power method, called the two-stage power method, applied to the matrix $A$ of part (a):

  1. Pick $\mathbf{x}^{(0)} \in \mathbb{R}^{n}$ satisfying $\left\|\mathbf{x}^{(0)}\right\|=1$. Let $0<\varepsilon \ll 1$. Set $k=0$ and $\mathbf{x}^{(1)}=A \mathbf{x}^{(0)}$.

  2. Calculate $\mathbf{x}^{(k+2)}=A \mathbf{x}^{(k+1)}$ and calculate $\alpha, \beta$ that minimise

  $$f(\alpha, \beta)=\left\|\mathbf{x}^{(k+2)}+\alpha \mathbf{x}^{(k+1)}+\beta \mathbf{x}^{(k)}\right\| .$$

  2. If $f(\alpha, \beta) \leqslant \varepsilon$, solve $\lambda^{2}+\alpha \lambda+\beta=0$ and let the roots be $\lambda_{1}$ and $\lambda_{2}$. They are accepted as eigenvalues of $A$, and the corresponding eigenvectors are estimated as $\mathbf{x}^{(k+1)}-\lambda_{2} \mathbf{x}^{(k)}$ and $\mathbf{x}^{(k+1)}-\lambda_{1} \mathbf{x}^{(k)}$.

  3. Otherwise, divide $\mathbf{x}^{(k+2)}$ and $\mathbf{x}^{(k+1)}$ by the current value of $\left\|\mathbf{x}^{(k+1)}\right\|$, increase $k$ by 1 and return to Step $\mathbf{1}$.

  Explain the justification behind Step 2 of the algorithm.

  (c) Let $n=3$, and suppose that, for a large value of $k$, the two-stage power method yields the vectors

  $$\mathbf{y}_{k}=\mathbf{x}^{(k)}=\left(\begin{array}{l} 1 \\ 1 \\ 1 \end{array}\right), \quad \mathbf{y}_{k+1}=A \mathbf{x}^{(k)}=\left(\begin{array}{c} 2 \\ 3 \\ 4 \end{array}\right), \quad \mathbf{y}_{k+2}=A^{2} \mathbf{x}^{(k)}=\left(\begin{array}{c} 2 \\ 4 \\ 6 \end{array}\right)$$

  Find two eigenvalues of $A$ and the corresponding eigenvectors.

  comment

• Paper 2, Section II, B

  Numerical Analysis

  Part II, 2016

  (a) The advection equation

  $$u_{t}=u_{x}, \quad 0 \leqslant x \leqslant 1, t \geqslant 0$$

  is discretised using an equidistant grid with stepsizes $\Delta x=h$ and $\Delta t=k$. The spatial derivatives are approximated with central differences and the resulting ODEs are approximated with the trapezoidal rule. Write down the relevant difference equation for determining $\left(u_{m}^{n+1}\right)$ from $\left(u_{m}^{n}\right)$. What is the name of this scheme? What is the local truncation error?

  The boundary condition is periodic, $u(0, t)=u(1, t)$. Explain briefly how to write the discretised scheme in the form $B \mathbf{u}^{n+1}=C \mathbf{u}^{n}$, where the matrices $B$ and $C$, to be identified, have a circulant form. Using matrix analysis, find the range of $\mu=\Delta t / \Delta x$ for which the scheme is stable. [Standard results may be used without proof if quoted carefully.]

  [Hint: An $n \times n$ circulant matrix has the form

  $$A=\left(\begin{array}{cccc} a_{0} & a_{1} & \cdots & a_{n-1} \\ a_{n-1} & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & a_{1} \\ a_{1} & \cdots & a_{n-1} & a_{0} \end{array}\right)$$

  All such matrices have the same set of eigenvectors $\mathbf{v}_{\ell}=\left(\omega^{j \ell}\right)_{j=0}^{n-1}, \ell=0,1, \ldots, n-1$, where $\omega=e^{2 \pi i / n}$, and the corresponding eigenvalues are $\lambda_{\ell}=\sum_{k=0}^{n-1} a_{k} \omega^{k \ell}$.]

  (b) Consider the advection equation on the unit square

  $$u_{t}=a u_{x}+b u_{y}, \quad 0 \leqslant x, y \leqslant 1, t \geqslant 0$$

  where $u$ satisfies doubly periodic boundary conditions, $u(0, y)=u(1, y), u(x, 0)=u(x, 1)$, and $a(x, y)$ and $b(x, y)$ are given doubly periodic functions. The system is discretised with the Crank-Nicolson scheme, with central differences for the space derivatives, using an equidistant grid with stepsizes $\Delta x=\Delta y=h$ and $\Delta t=k$. Write down the relevant difference equation, and show how to write the scheme in the form

  $$\mathbf{u}^{n+1}=\left(I-\frac{1}{4} \mu A\right)^{-1}\left(I+\frac{1}{4} \mu A\right) \mathbf{u}^{n}$$

  where the matrix $A$ should be identified. Describe how (*) can be approximated by Strang splitting, and explain the advantages of doing so.

  [Hint: Inversion of the matrix $B$ in part (a) has a similar computational cost to that of a tridiagonal matrix.]

  comment

• Paper 1, Section II, B

  Numerical Analysis

  Part II, 2016

  (a) Consider the periodic function

  $$f(x)=5+2 \cos \left(2 \pi x-\frac{\pi}{2}\right)+3 \cos (4 \pi x)$$

  on the interval $[0,1]$. The $N$-point discrete Fourier transform of $f$ is defined by

  $$F_{N}(n)=\frac{1}{N} \sum_{k=0}^{N-1} f_{k} \omega_{N}^{-n k}, \quad n=0,1, \ldots, N-1$$

  where $\omega_{N}=e^{2 \pi i / N}$ and $f_{k}=f(k / N)$. Compute $F_{4}(n), n=0, \ldots, 3$, using the Fast Fourier Transform (FFT). Compare your result with what you get by computing $F_{4}(n)$ directly from $(*)$. Carefully explain all your computations.

  (b) Now let $f$ be any analytic function on $\mathbb{R}$ that is periodic with period 1 . The continuous Fourier transform of $f$ is defined by

  $$\hat{f}_{n}=\int_{0}^{1} f(\tau) e^{-2 \pi i n \tau} d \tau, \quad n \in \mathbb{Z}$$

  Use integration by parts to show that the Fourier coefficients $\hat{f}_{n}$ decay spectrally.

  Explain what it means for the discrete Fourier transform of $f$ to approximate the continuous Fourier transform with spectral accuracy. Prove that it does so.

  What can you say about the behaviour of $F_{N}(N-n)$ as $N \rightarrow \infty$ for fixed $n$ ?

  comment

• Paper 3, Section II, B

  Numerical Analysis

  Part II, 2016

  (a) Define the Jacobi and Gauss-Seidel iteration schemes for solving a linear system of the form $A \mathbf{u}=\mathbf{b}$, where $\mathbf{u}, \mathbf{b} \in \mathbb{R}^{M}$ and $A \in \mathbb{R}^{M \times M}$, giving formulae for the corresponding iteration matrices $H_{J}$ and $H_{G S}$ in terms of the usual decomposition $A=L_{0}+D+U_{0}$.

  Show that both iteration schemes converge when $A$ results from discretization of Poisson's equation on a square with the five-point formula, that is when

  $$A=\left[\begin{array}{rrrrr} S & I & & & \\ I & S & I & & \\ & \ddots & \ddots & \ddots & \\ & & I & S & I \\ & & & I & S \end{array}\right], \quad S=\left[\begin{array}{rrrrr} -4 & 1 & & & \\ 1 & -4 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -4 & 1 \\ & & & 1 & -4 \end{array}\right] \in \mathbb{R}^{m \times m}$$

  and $M=m^{2}$. [You may state the Householder-John theorem without proof.]

  (b) For the matrix $A$ given in $(*)$ :

  (i) Calculate the eigenvalues of $H_{J}$ and deduce its spectral radius $\rho\left(H_{J}\right)$.

  (ii) Show that each eigenvector $\mathbf{q}$ of $H_{G S}$ is related to an eigenvector $\mathbf{p}$ of $H_{J}$ by a transformation of the form $q_{i, j}=\alpha^{i+j} p_{i, j}$ for a suitable value of $\alpha$.

  (iii) Deduce that $\rho\left(H_{G S}\right)=\rho^{2}\left(H_{J}\right)$. What is the significance of this result for the two iteration schemes?

  [Hint: You may assume that the eigenvalues of the matrix $A$ in $(*)$ are

  $$\lambda_{k, \ell}=-4\left(\sin ^{2} \frac{x}{2}+\sin ^{2} \frac{y}{2}\right), \quad \text { where } x=\frac{k \pi h}{m+1}, y=\frac{\ell \pi h}{m+1}, \quad k, \ell=1, \ldots, m,$$

  with corresponding eigenvectors $\left.\mathbf{v}=\left(v_{i, j}\right), \quad v_{i, j}=\sin i x \sin j y, \quad i, j=1, \ldots, m . \quad\right]$

  comment

• Paper 2, Section II, A

  Numerical Analysis

  Part II, 2017

  The Poisson equation $\nabla^{2} u=f$ in the unit square $\Omega=[0,1] \times[0,1]$, equipped with the zero Dirichlet boundary conditions on $\partial \Omega$, is discretized with the nine-point formula:

  $$\begin{aligned} \Gamma_{9}\left[u_{i, j}\right]:=&-\frac{10}{3} u_{i, j}+\frac{2}{3}\left(u_{i+1, j}+u_{i-1, j}+u_{i, j+1}+u_{i, j-1}\right) \\ &+\frac{1}{6}\left(u_{i+1, j+1}+u_{i+1, j-1}+u_{i-1, j+1}+u_{i-1, j-1}\right)=h^{2} f_{i, j} \end{aligned}$$

  where $1 \leqslant i, j \leqslant m, u_{i, j} \approx u(i h, j h)$, and $(i h, j h)$ are the grid points with $h=\frac{1}{m+1}$.

  (i) Find the order of the local truncation error $\eta_{i, j}$ of the approximation.

  (ii) Prove that the order of the truncation error is smaller if $f$ satisfies the Laplace equation $\nabla^{2} f=0$.

  (iii) Show that the modified nine-point scheme

  $$\begin{aligned} \Gamma_{9}\left[u_{i, j}\right] &=h^{2} f_{i, j}+\frac{1}{12} h^{2} \Gamma_{5}\left[f_{i, j}\right] \\ &:=h^{2} f_{i, j}+\frac{1}{12} h^{2}\left(f_{i+1, j}+f_{i-1, j}+f_{i, j+1}+f_{i, j-1}-4 f_{i, j}\right) \end{aligned}$$

  has a truncation error of the same order as in part (ii).

  (iv) Let $\left(u_{i, j}\right)_{i, j=1}^{m}$ be a solution to the $m^{2} \times m^{2}$ system of linear equations $A \mathbf{u}=\mathbf{b}$ arising from the modified nine-point scheme in part (iii). Further, let $u(x, y)$ be the exact solution and let $e_{i, j}:=u_{i, j}-u(i h, j h)$ be the error of approximation at grid points. Prove that there exists a constant $c$ such that

  $$\|\mathbf{e}\|_{2}:=\left[\sum_{i, j=1}^{m}\left|e_{i, j}\right|^{2}\right]^{1 / 2}<c h^{3}, \quad h \rightarrow 0$$

  [Hint: The nine-point discretization of $\nabla^{2} u$ can be written as

  $$\Gamma_{9}[u]=\left(\Gamma_{5}+\frac{1}{6} \Delta_{x}^{2} \Delta_{y}^{2}\right) u=\left(\Delta_{x}^{2}+\Delta_{y}^{2}+\frac{1}{6} \Delta_{x}^{2} \Delta_{y}^{2}\right) u$$

  where $\Gamma_{5}[u]=\left(\Delta_{x}^{2}+\Delta_{y}^{2}\right) u$ is the five-point discretization and

  $$\left.\begin{array}{l} \Delta_{x}^{2} u(x, y):=u(x-h, y)-2 u(x, y)+u(x+h, y) \\ \Delta_{y}^{2} u(x, y):=u(x, y-h)-2 u(x, y)+u(x, y+h) \end{array}\right]$$

  [Hint: The matrix A of the nine-point scheme is symmetric, with the eigenvalues

  $$\lambda_{k, \ell}=-4 \sin ^{2} \frac{k \pi h}{2}-4 \sin ^{2} \frac{\ell \pi h}{2}+\frac{8}{3} \sin ^{2} \frac{k \pi h}{2} \sin ^{2} \frac{\ell \pi h}{2}, \quad 1 \leqslant k, \ell \leqslant m$$

  comment

• Paper 1, Section II, A

  Numerical Analysis

  Part II, 2017

  State the Householder-John theorem and explain how it can be used in designing iterative methods for solving a system of linear equations $A \mathbf{x}=\mathbf{b}$. [You may quote other relevant theorems if needed.]

  Consider the following iterative scheme for solving $A \mathbf{x}=\mathbf{b}$. Let $A=L+D+U$, where $D$ is the diagonal part of $A$, and $L$ and $U$ are the strictly lower and upper triangular parts of $A$, respectively. Then, with some starting vector $\mathbf{x}^{(0)}$, the scheme is as follows:

  $$(D+\omega L) \mathbf{x}^{(k+1)}=[(1-\omega) D-\omega U] \mathbf{x}^{(k)}+\omega \mathbf{b} .$$

  Prove that if $A$ is a symmetric positive definite matrix and $\omega \in(0,2)$, then, for any $\mathbf{x}^{(0)}$, the above iteration converges to the solution of the system $A \mathbf{x}=\mathbf{b}$.

  Which method corresponds to the case $\omega=1 ?$

  comment

• Paper 3, Section II, A

  Numerical Analysis

  Part II, 2017

  Let $A$ be a real symmetric $n \times n$ matrix with real and distinct eigenvalues $0=\lambda_{1}<\cdots<\lambda_{n-1}=1<\lambda_{n}$ and a corresponding orthogonal basis of normalized real eigenvectors $\left(\mathbf{w}_{i}\right)_{i=1}^{n}$.

  To estimate the eigenvector $\mathbf{w}_{n}$ of $A$ whose eigenvalue is $\lambda_{n}$, the power method with shifts is employed which has the following form:

  $$\mathbf{y}=\left(A-s_{k} I\right) \mathbf{x}^{(k)}, \quad \mathbf{x}^{(k+1)}=\mathbf{y} /\|\mathbf{y}\|, \quad s_{k} \in \mathbb{R}, \quad k=0,1,2, \ldots$$

  Three versions of this method are considered:

  (i) no shift: $s_{k} \equiv 0$;

  (ii) single shift: $s_{k} \equiv \frac{1}{2}$;

  (iii) double shift: $s_{2 \ell} \equiv s_{0}=\frac{1}{4}(2+\sqrt{2}), s_{2 \ell+1} \equiv s_{1}=\frac{1}{4}(2-\sqrt{2})$.

  Assume that $\lambda_{n}=1+\epsilon$, where $\epsilon>0$ is very small, so that the terms $\mathcal{O}\left(\epsilon^{2}\right)$ are negligible, and that $\mathbf{x}^{(0)}$ contains substantial components of all the eigenvectors.

  By considering the approximation after $2 m$ iterations in the form

  $$\mathbf{x}^{(2 m)}=\pm \mathbf{w}_{n}+\mathcal{O}\left(\rho^{2 m}\right) \quad(m \rightarrow \infty)$$

  find $\rho$ as a function of $\epsilon$ for each of the three versions of the method.

  Compare the convergence rates of the three versions of the method, with reference to the number of iterations needed to achieve a prescribed accuracy.

  comment

• Paper 4, Section II, A

  Numerical Analysis

  Part II, 2017

  (a) The diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0 \leqslant x \leqslant 1, \quad 0 \leqslant t \leqslant T$$

  is approximated by the Crank-Nicolson scheme

  $$u_{m}^{n+1}-\frac{1}{2} \mu\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right)=u_{m}^{n}+\frac{1}{2} \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)$$

  with $m=1, \ldots, M$. Here $\mu=k / h^{2}, k=\Delta t, h=\Delta x=\frac{1}{M+1}$, and $u_{m}^{n}$ is an approximation to $u(m h, n k)$. Assuming that $u(0, t)=u(1, t) \stackrel{M+1}{=} 0$, show that the above scheme can be written in the form

  $$B \mathbf{u}^{n+1}=C \mathbf{u}^{n}, \quad 0 \leqslant n \leqslant T / k-1$$

  where $\mathbf{u}^{n}=\left[u_{1}^{n}, \ldots, u_{M}^{n}\right]^{T}$ and the real matrices $B$ and $C$ should be found. Using matrix analysis, find the range of $\mu>0$ for which the scheme is stable.

  [Hint: All Toeplitz symmetric tridiagonal (TST) matrices have the same set of orthogonal eigenvectors, and a TST matrix with the elements $a_{i, i}=a$ and $a_{i, i \pm 1}=b$ has the eigenvalues $\lambda_{k}=a+2 b \cos \frac{\pi k}{M+1}$. ]

  (b) The wave equation

  $$\frac{\partial^{2} u}{\partial t^{2}}=\frac{\partial^{2} u}{\partial x^{2}}, \quad x \in \mathbb{R}, \quad t \geqslant 0$$

  with given initial conditions for $u$ and $\partial u / \partial t$, is approximated by the scheme

  $$u_{m}^{n+1}-2 u_{m}^{n}+u_{m}^{n-1}=\mu\left(u_{m+1}^{n}-2 u_{m}^{n}+u_{m-1}^{n}\right),$$

  with the Courant number now $\mu=k^{2} / h^{2}$. Applying the Fourier technique, find the range of $\mu>0$ for which the method is stable.

  comment

• Paper 4, Section II, E

  Numerical Analysis

  Part II, 2018

  The inverse discrete Fourier transform $\mathcal{F}_{n}^{-1}: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ is given by the formula

  $$\boldsymbol{x}=\mathcal{F}_{n}^{-1} \boldsymbol{y}, \quad \text { where } \quad x_{\ell}=\sum_{j=0}^{n-1} \omega_{n}^{j \ell} y_{j}, \quad \ell=0, \ldots, n-1$$

  Here, $\omega_{n}=\exp \frac{2 \pi i}{n}$ is the primitive root of unity of degree $n$ and $n=2^{p}, p=1,2, \ldots$

  (a) Show how to assemble $\boldsymbol{x}=\mathcal{F}_{2 m}^{-1} \boldsymbol{y}$ in a small number of operations if the Fourier transforms of the even and odd parts of $\boldsymbol{y}$,

  $$\boldsymbol{x}^{(\mathrm{E})}=\mathcal{F}_{m}^{-1} \boldsymbol{y}^{(\mathrm{E})}, \quad \boldsymbol{x}^{(\mathrm{O})}=\mathcal{F}_{m}^{-1} \boldsymbol{y}^{(\mathrm{O})}$$

  are already known.

  (b) Describe the Fast Fourier Transform (FFT) method for evaluating $\boldsymbol{x}$, and draw a diagram to illustrate the method for $n=8$.

  (c) Find the cost of the FFT method for $n=2^{p}$ (only multiplications count).

  (d) For $n=4$ use the FFT method to find $\boldsymbol{x}=\mathcal{F}_{n}^{-1} \boldsymbol{y}$ when: (i) $\boldsymbol{y}=(1,-1,1,-1)$, (ii) $\boldsymbol{y}=(1,1,-1,-1)$.

  comment

• Paper 2, Section II, E

  Numerical Analysis

  Part II, 2018

  The Poisson equation $\frac{d^{2} u}{d x^{2}}=f$ in the unit interval $[0,1]$, with $u(0)=u(1)=0$, is discretised with the formula

  $$u_{i-1}-2 u_{i}+u_{i+1}=h^{2} f_{i}, \quad 1 \leqslant i \leqslant n,$$

  where $u_{0}=u_{n+1}=0, h=(n+1)^{-1}$, the grid points are at $x=i h$ and $u_{i} \approx u(i h)$.

  (a) Write the above system of equations in the vector form $A \boldsymbol{u}=\boldsymbol{b}$ and describe the relaxed Jacobi method with relaxation parameter $\omega$ for solving this linear system.

  (b) For $\boldsymbol{x}^{*}$ and $\boldsymbol{x}^{(\nu)}$ being the exact and the iterated solution, respectively, let $\boldsymbol{e}^{(\nu)}:=\boldsymbol{x}^{(\nu)}-\boldsymbol{x}^{*}$ be the error and $H_{\omega}$ be the iteration matrix, so that

  $$\boldsymbol{e}^{(\nu+1)}=H_{\omega} \boldsymbol{e}^{(\nu)}$$

  Express $H_{\omega}$ in terms of the matrix $A$ and the relaxation parameter $\omega$. Using the fact that for any $n \times n$ Toeplitz symmetric tridiagonal matrix, the eigenvectors $\boldsymbol{v}_{k}(k=1, \ldots, n)$ have the form:

  $$\boldsymbol{v}_{k}=(\sin i k x)_{i=1}^{n}, \quad x=\frac{\pi}{n+1}$$

  find the eigenvalues $\lambda_{k}(A)$ of $A$. Hence deduce the eigenvalues $\lambda_{k}(\omega)$ of $H_{\omega}$.

  (c) For $A$ as above, let

  $$\boldsymbol{e}^{(\nu)}=\sum_{k=1}^{n} a_{k}^{(\nu)} \boldsymbol{v}_{k}$$

  be the expansion of the error with respect to the eigenvectors $\left(\boldsymbol{v}_{k}\right)$ of $H_{\omega}$.

  Find the range of the parameter $\omega$ which provides convergence of the method for any $n$, and prove that, for any such $\omega$, the rate of convergence $e^{(\nu)} \rightarrow 0$ is not faster than $\left(1-c / n^{2}\right)^{\nu}$ when $n$ is large.

  (d) Show that, for an appropriate range of $\omega$, the high frequency components $a_{k}^{(\nu)}$ $\left(\frac{n+1}{2} \leqslant k \leqslant n\right)$ of the error $e^{(\nu)}$ tend to zero much faster than the rate obtained in part (c). Determine the optimal parameter $\omega_{*}$ which provides the largest supression of the high frequency components per iteration, and find the corresponding attenuation factor $\mu_{*}$ assuming $n$ is large. That is, find the least $\mu_{\omega}$ such that $\left|a_{k}^{(\nu+1)}\right| \leqslant \mu_{\omega}\left|a_{k}^{(\nu)}\right|$ for $\frac{n+1}{2} \leqslant k \leqslant n$.

  comment

• Paper 1, Section II, E

  Numerical Analysis

  Part II, 2018

  (a) Suppose that $A$ is a real $n \times n$ matrix, and $\boldsymbol{w} \in \mathbb{R}^{n}$ and $\lambda_{1} \in \mathbb{R}$ are given so that $A \boldsymbol{w}=\lambda_{1} \boldsymbol{w}$. Further, let $S$ be a non-singular matrix such that $S \boldsymbol{w}=c e^{(1)}$, where $e^{(1)}$ is the first coordinate vector and $c \neq 0$.

  Let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are $\lambda_{1}$ together with the eigenvalues of the bottom right $(n-1) \times(n-1)$ submatrix of $\widehat{A}$.

  Explain briefly how, given a vector $\boldsymbol{w}$, an orthogonal matrix $S$ such that $S \boldsymbol{w}=c e^{(1)}$ can be constructed.

  (b) Suppose that $A$ is a real $n \times n$ matrix, and two linearly independent vectors $\boldsymbol{v}, \boldsymbol{w} \in \mathbb{R}^{n}$ are given such that the linear subspace $L\{\boldsymbol{v}, \boldsymbol{w}\}$ spanned by $\boldsymbol{v}$ and $\boldsymbol{w}$ is invariant under the action of $A$, i.e.,

  $$x \in L\{\boldsymbol{v}, \boldsymbol{w}\} \quad \Rightarrow \quad A x \in L\{\boldsymbol{v}, \boldsymbol{w}\}$$

  Denote by $V$ an $n \times 2$ matrix whose two columns are the vectors $\boldsymbol{v}$ and $\boldsymbol{w}$, and let $S$ be a non-singular matrix such that $R=S V$ is upper triangular:

  $$S V=S \times\left[\begin{array}{cc} v_{1} & w_{1} \\ v_{2} & w_{2} \\ v_{3} & w_{3} \\ : & : \\ v_{n} & w_{n} \end{array}\right]=\left[\begin{array}{cc} r_{11} & r_{12} \\ 0 & r_{22} \\ 0 & 0 \\ : & : \\ 0 & 0 \end{array}\right]$$

  Again, let $\widehat{A}=S A S^{-1}$. Prove that the eigenvalues of $A$ are the eigenvalues of the top left $2 \times 2$ submatrix of $\widehat{A}$together with the eigenvalues of the bottom right $(n-2) \times(n-2)$ submatrix of $\widehat{A}$.

  Explain briefly how, for given vectors $\boldsymbol{v}, \boldsymbol{w}$, an orthogonal matrix $S$ which satisfies (*) can be constructed.

  comment

• Paper 3, Section II, E

  Numerical Analysis

  Part II, 2018

  The diffusion equation for $u(x, t)$ :

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad x \in \mathbb{R}, \quad t \geqslant 0$$

  is solved numerically by the difference scheme

  $$u_{m}^{n+1}=u_{m}^{n}+\frac{3}{2} \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)-\frac{1}{2} \mu\left(u_{m-1}^{n-1}-2 u_{m}^{n-1}+u_{m+1}^{n-1}\right) .$$

  Here $\mu=\frac{k}{h^{2}}$ is the Courant number, with $k=\Delta t, h=\Delta x$, and $u_{m}^{n} \approx u(m h, n k)$.

  (a) Prove that, as $k \rightarrow 0$ with constant $\mu$, the local error of the method is $\mathcal{O}\left(k^{2}\right)$.

  (b) Applying the Fourier stability analysis, show that the method is stable if and only if $\mu \leqslant \frac{1}{4}$. [Hint: If a polynomial $p(x)=x^{2}-2 \alpha x+\beta$ has real roots, then those roots lie in $[a, b]$ if and only if $p(a) p(b) \geqslant 0$ and $\alpha \in[a, b]$.]

  (c) Prove that, for the same equation, the leapfrog scheme

  $$u_{m}^{n+1}=u_{m}^{n-1}+2 \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)$$

  is unstable for any choice of $\mu>0$.

  comment

• Paper 4, Section II, C

  Numerical Analysis

  Part II, 2019

  For a 2-periodic analytic function $f$, its Fourier expansion is given by the formula

  $$f(x)=\sum_{n=-\infty}^{\infty} \widehat{f}_{n} e^{i \pi n x}, \quad \widehat{f}_{n}=\frac{1}{2} \int_{-1}^{1} f(t) e^{-i \pi n t} d t$$

  (a) Consider the two-point boundary value problem

  $$-\frac{1}{\pi^{2}}(1+2 \cos \pi x) u^{\prime \prime}+u=1+\sum_{n=1}^{\infty} \frac{2}{n^{2}+1} \cos \pi n x, \quad-1 \leqslant x \leqslant 1,$$

  with periodic boundary conditions $u(-1)=u(1)$. Construct explicitly the infinite dimensional linear algebraic system that arises from the application of the Fourier spectral method to the above equation, and explain how to truncate the system to a finitedimensional one.

  (b) A rectangle rule is applied to computing the integral of a 2-periodic analytic function $h$ :

  $$\int_{-1}^{1} h(t) d t \approx \frac{2}{N} \sum_{k=-N / 2+1}^{N / 2} h\left(\frac{2 k}{N}\right)$$

  Find an expression for the error $e_{N}(h):=\mathrm{LHS}-\mathrm{RHS}$ of $(*)$, in terms of $\widehat{h}_{n}$, and show that $e_{N}(h)$ has a spectral rate of decay as $N \rightarrow \infty$. [In the last part, you may quote a relevant theorem about $\widehat{h}_{n}$.]

  comment

• Paper 2, Section II, C

  Numerical Analysis

  Part II, 2019

  The Poisson equation on the unit square, equipped with zero boundary conditions, is discretized with the 9-point scheme:

  $$\begin{aligned} &-\frac{10}{3} u_{i, j}+\frac{2}{3}\left(u_{i+1, j}+u_{i-1, j}+u_{i, j+1}+u_{i, j-1}\right) \\ &+\frac{1}{6}\left(u_{i+1, j+1}+u_{i+1, j-1}+u_{i-1, j+1}+u_{i-1, j-1}\right)=h^{2} f_{i, j} \end{aligned}$$

  where $1 \leqslant i, j \leqslant m, u_{i, j} \approx u(i h, j h)$, and $(i h, j h)$ are the grid points with $h=\frac{1}{m+1}$. We also assume that $u_{0, j}=u_{i, 0}=u_{m+1, j}=u_{i, m+1}=0$.

  (a) Prove that all $m \times m$ tridiagonal symmetric Toeplitz (TST-) matrices

  $$H=[\beta, \alpha, \beta]:=\left[\begin{array}{cccc} \alpha & \beta & & \\ \beta & \alpha & \ddots & \\ & \ddots & \ddots & \beta \\ & & \beta & \alpha \end{array}\right] \in \mathbb{R}^{m \times m}$$

  share the same eigenvectors $\boldsymbol{q}_{k}$ with the components $(\sin j k \pi h)_{j=1}^{m}$ for $k=1, \ldots, m$. Find expressions for the corresponding eigenvalues $\lambda_{k}$ for $k=1, \ldots, m$. Deduce that $H=Q D Q^{-1}$, where $D=\operatorname{diag}\left\{\lambda_{k}\right\}$ and $Q$ is the matrix $(\sin i j \pi h)_{i, j=1}^{m} .$

  (b) Show that, by arranging the grid points ( $i h, j h)$ into a one-dimensional array by columns, the 9 -points scheme results in the following system of linear equations of the form

  $$A \boldsymbol{u}=\boldsymbol{b} \Leftrightarrow\left[\begin{array}{cccc} B & C & & \\ C & B & \ddots & \\ & \ddots & \ddots & C \\ & & C & B \end{array}\right]\left[\begin{array}{c} \boldsymbol{u}_{1} \\ \boldsymbol{u}_{2} \\ \vdots \\ \boldsymbol{u}_{m} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{b}_{1} \\ \boldsymbol{b}_{2} \\ \vdots \\ \boldsymbol{b}_{m} \end{array}\right]$$

  where $A \in \mathbb{R}^{m^{2} \times m^{2}}$, the vectors $\boldsymbol{u}_{k}, \boldsymbol{b}_{k} \in \mathbb{R}^{m}$ are portions of $\boldsymbol{u}, \boldsymbol{b} \in \mathbb{R}^{m^{2}}$, respectively, and $B, C$ are $m \times m$ TST-matrices whose elements you should determine.

  (c) Using $\boldsymbol{v}_{k}=Q^{-1} \boldsymbol{u}_{k}, \boldsymbol{c}_{k}=Q^{-1} \boldsymbol{b}_{k}$, show that (2) is equivalent to

  $$\left[\begin{array}{cccc} D & E & & \\ E & D & \ddots & \\ & \ddots & \ddots & E \\ & & E & D \end{array}\right]\left[\begin{array}{c} \boldsymbol{v}_{1} \\ \boldsymbol{v}_{2} \\ \vdots \\ \boldsymbol{v}_{m} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{c}_{1} \\ \boldsymbol{c}_{2} \\ \vdots \\ \boldsymbol{c}_{m} \end{array}\right]$$

  where $D$ and $E$ are diagonal matrices.

  (d) Show that, by appropriate reordering of the grid, the system (3) is reduced to $m$ uncoupled $m \times m$ systems of the form

  $$\Lambda_{k} \widehat{v}_{k}=\widehat{c}_{k}, \quad k=1, \ldots, m$$

  Determine the elements of the matrices $\Lambda_{k}$.

  comment

• Paper 3, Section II, 40C

  Numerical Analysis

  Part II, 2019

  The diffusion equation

  $$u_{t}=u_{x x}, \quad 0 \leqslant x \leqslant 1, \quad t \geqslant 0,$$

  with the initial condition $u(x, 0)=\phi(x), 0 \leqslant x \leqslant 1$, and boundary conditions $u(0, t)=$ $u(1, t)=0$, is discretised by $u_{m}^{n} \approx u(m h, n k)$ with $k=\Delta t, h=\Delta x=1 /(1+M)$. The Courant number is given by $\mu=k / h^{2}$.

  (a) The system is solved numerically by the method

  $$u_{m}^{n+1}=u_{m}^{n}+\mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right), \quad m=1,2, \ldots, M, \quad n \geqslant 0 .$$

  Prove directly that $\mu \leqslant 1 / 2$ implies convergence.

  (b) Now consider the method

  $$a u_{m}^{n+1}-\frac{1}{4}(\mu-c)\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right)=a u_{m}^{n}+\frac{1}{4}(\mu+c)\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right)$$

  where $a$ and $c$ are real constants. Using an eigenvalue analysis and carefully justifying each step, determine conditions on $\mu, a$ and $c$ for this method to be stable.

  [You may use the notation $[\beta, \alpha, \beta]$ for the tridiagonal matrix with $\alpha$ along the diagonal, and $\beta$ along the sub-and super-diagonals and use without proof any relevant theorems about such matrices.]

  comment

• Paper 1, Section II, C

  Numerical Analysis

  Part II, 2019

  (a) Describe the Jacobi method for solving a system of linear equations $A \boldsymbol{x}=\boldsymbol{b}$ as a particular case of splitting, and state the criterion for its convergence in terms of the iteration matrix.

  (b) For the case when

  $$A=\left[\begin{array}{lll} 1 & \alpha & \alpha \\ \alpha & 1 & \alpha \\ \alpha & \alpha & 1 \end{array}\right]$$

  find the exact range of the parameter $\alpha$ for which the Jacobi method converges.

  (c) State the Householder-John theorem and deduce that the Jacobi method converges if $A$ is a symmetric positive-definite tridiagonal matrix.

  comment

• Paper 1, Section II, E

  Numerical Analysis

  Part II, 2020

  Let $A \in \mathbb{R}^{n \times n}$ be a real symmetric matrix with distinct eigenvalues $\lambda_{1}<\lambda_{2}<\cdots<$ $\lambda_{n}$ and a corresponding orthonormal basis of real eigenvectors $\left\{\boldsymbol{w}_{i}\right\}_{i=1}^{n}$. Given a unit norm vector $\boldsymbol{x}^{(0)} \in \mathbb{R}^{n}$, and a set of parameters $s_{k} \in \mathbb{R}$, consider the inverse iteration algorithm

  $$\left(A-s_{k} I\right) \boldsymbol{y}=\boldsymbol{x}^{(k)}, \quad \boldsymbol{x}^{(k+1)}=\boldsymbol{y} /\|\boldsymbol{y}\|, \quad k \geqslant 0 .$$

  (a) Let $s_{k}=s=$ const for all $k$. Assuming that $\boldsymbol{x}^{(0)}=\sum_{i=1}^{n} c_{i} \boldsymbol{w}_{i}$ with all $c_{i} \neq 0$, prove that

  $$s<\lambda_{1} \Rightarrow \boldsymbol{x}^{(k)} \rightarrow \boldsymbol{w}_{1} \quad \text { or } \quad \boldsymbol{x}^{(k)} \rightarrow-\boldsymbol{w}_{1} \quad(k \rightarrow \infty) .$$

  Explain briefly what happens to $\boldsymbol{x}^{(k)}$ when $\lambda_{m}<s<\lambda_{m+1}$ for some $m \in\{1,2, \ldots, n-1\}$, and when $\lambda_{n}<s$.

  (b) Let $s_{k}=\left(A \boldsymbol{x}^{(k)}, \boldsymbol{x}^{(k)}\right)$ for $k \geqslant 0$. Assuming that, for some $k$, some $a_{i} \in \mathbb{R}$ and for a small $\epsilon$,

  $$\boldsymbol{x}^{(k)}=c^{-1}\left(\boldsymbol{w}_{1}+\epsilon \sum_{i \geqslant 2} a_{i} \boldsymbol{w}_{i}\right)$$

  where $c$ is the appropriate normalising constant. Show that $s_{k}=\lambda_{1}-K \epsilon^{2}+\mathcal{O}\left(\epsilon^{4}\right)$ and determine the value of $K$. Hence show that

  $$\boldsymbol{x}^{(k+1)}=c_{1}^{-1}\left(\boldsymbol{w}_{1}+\epsilon^{3} \sum_{i \geqslant 2} b_{i} \boldsymbol{w}_{i}+\mathcal{O}\left(\epsilon^{5}\right)\right)$$

  where $c_{1}$ is the appropriate normalising constant, and find expressions for $b_{i}$.

  comment

• Paper 2, Section II, 40E

  Numerical Analysis

  Part II, 2020

  (a) For $A \in \mathbb{R}^{n \times n}$ and nonzero $\boldsymbol{v} \in \mathbb{R}^{n}$, define the $m$-th Krylov subspace $K_{m}(A, \boldsymbol{v})$ of $\mathbb{R}^{n}$. Prove that if $A$ has $n$ linearly independent eigenvectors with at most $s$ distinct eigenvalues, then

  $$\operatorname{dim} K_{m}(A, \boldsymbol{v}) \leqslant s \quad \forall m$$

  (b) Define the term residual in the conjugate gradient (CG) method for solving a system $A \boldsymbol{x}=\boldsymbol{b}$ with a symmetric positive definite $A$. State two properties of the method regarding residuals and their connection to certain Krylov subspaces, and hence show that, for any right-hand side $\boldsymbol{b}$, the method finds the exact solution after at most $s$ iterations, where $s$ is the number of distinct eigenvalues of $A$.

  (c) The preconditioned CG-method $P A P^{T} \widehat{\boldsymbol{x}}=P \boldsymbol{b}$ is applied for solving $A \boldsymbol{x}=\boldsymbol{b}$, with

  

  Prove that the method finds the exact solution after two iterations at most.

  (d) Prove that, for any symmetric positive definite $A$, we can find a preconditioner $P$ such that the preconditioned CG-method for solving $A \boldsymbol{x}=\boldsymbol{b}$ would require only one step. Explain why this preconditioning is of hardly any use.

  comment

• Paper 3, Section II, 40E

  Numerical Analysis

  Part II, 2020

  (a) Give the definition of a normal matrix. Prove that if $A$ is normal, then the (Euclidean) matrix $\ell_{2}$-norm of $A$ is equal to its spectral radius, i.e., $\|A\|_{2}=\rho(A)$.

  (b) The advection equation

  $$u_{t}=u_{x}, \quad 0 \leqslant x \leqslant 1, \quad 0 \leqslant t<\infty$$

  is discretized by the Crank-Nicolson scheme

  $$u_{m}^{n+1}-u_{m}^{n}=\frac{1}{4} \mu\left(u_{m+1}^{n+1}-u_{m-1}^{n+1}\right)+\frac{1}{4} \mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right), \quad m=1,2, \ldots, M, \quad n \in \mathbb{Z}_{+}$$

  Here, $\mu=\frac{k}{h}$ is the Courant number, with $k=\Delta t, h=\Delta x=\frac{1}{M+1}$, and $u_{m}^{n}$ is an approximation to $u(m h, n k)$.

  Using the eigenvalue analysis and carefully justifying each step, determine conditions on $\mu>0$ for which the method is stable. [Hint: All M $\times M$ Toeplitz anti-symmetric tridiagonal (TAT) matrices have the same set of orthogonal eigenvectors, and a TAT matrix with the elements $a_{j, j}=a$ and $a_{j, j+1}=-a_{j, j-1}=b$ has the eigenvalues $\lambda_{k}=a+2 \mathrm{i} b \cos \frac{\pi k}{M+1}$ where $\mathrm{i}=\sqrt{-1}$. ]

  (c) Consider the same advection equation for the Cauchy problem $(x \in \mathbb{R}, 0 \leqslant t \leqslant$ $T)$. Now it is discretized by the two-step leapfrog scheme

  $$u_{m}^{n+1}=\mu\left(u_{m+1}^{n}-u_{m-1}^{n}\right)+u_{m}^{n-1} .$$

  Applying the Fourier technique, find the range of $\mu>0$ for which the method is stable.

  comment

• Paper 4 , Section II, 40E

  Numerical Analysis

  Part II, 2020

  (a) For a function $f=f(x, y)$ which is real analytic in $\mathbb{R}^{2}$ and 2-periodic in each variable, its Fourier expansion is given by the formula

  $$f(x, y)=\sum_{m, n \in \mathbb{Z}} \widehat{f}_{m, n} e^{i \pi m x+i \pi n y}, \quad \widehat{f}_{m, n}=\frac{1}{4} \int_{-1}^{1} \int_{-1}^{1} f(t, \theta) e^{-i \pi m t-i \pi n \theta} d t d \theta$$

  Derive expressions for the Fourier coefficients of partial derivatives $f_{x}, f_{y}$ and those of the product $h(x, y)=f(x, y) g(x, y)$ in terms of $\widehat{f}_{m, n}$ and $\widehat{g}_{m, n}$.

  (b) Let $u(x, y)$ be the 2-periodic solution in $\mathbb{R}^{2}$ of the general second-order elliptic PDE

  $$\left(a u_{x}\right)_{x}+\left(a u_{y}\right)_{y}=f$$

  where $a$ and $f$ are both real analytic and 2-periodic, and $a(x, y)>0$. We impose the normalisation condition $\int_{-1}^{1} \int_{-1}^{1} u d x d y=0$ and note from the PDE $\int_{-1}^{1} \int_{-1}^{1} f d x d y=0$.

  Construct explicitly the infinite-dimensional linear algebraic system that arises from the application of the Fourier spectral method to the above equation, and explain how to truncate this system to a finite-dimensional one.

  (c) Specify the truncated system for the unknowns $\left\{\widehat{u}_{m, n}\right\}$ for the case

  $$a(x, y)=5+2 \cos \pi x+2 \cos \pi y$$

  and prove that, for any ordering of the Fourier coefficients $\left\{\widehat{u}_{m, n}\right\}$ into one-dimensional array, the resulting system is symmetric and positive definite. [You may use the Gershgorin theorem without proof.]

  comment

• Paper 1, Section II, E

  Numerical Analysis

  Part II, 2021

  Let $A \in \mathbb{R}^{n \times n}$ with $n>2$ and define $\operatorname{Spec}(A)=\{\lambda \in \mathbb{C} \mid A-\lambda I$ is not invertible $\}$.

  The QR algorithm for computing $\operatorname{Spec}(A)$ is defined as follows. Set $A_{0}=A$. For $k=0,1, \ldots$ compute the $\mathrm{QR}$ factorization $A_{k}=Q_{k} R_{k}$ and set $A_{k+1}=R_{k} Q_{k}$. (Here $Q_{k}$ is an $n \times n$ orthogonal matrix and $R_{k}$ is an $n \times n$ upper triangular matrix.)

  (a) Show that $A_{k+1}$ is related to the original matrix $A$ by the similarity transformation $A_{k+1}=\bar{Q}_{k}^{T} A \bar{Q}_{k}$, where $\bar{Q}_{k}=Q_{0} Q_{1} \cdots Q_{k}$ is orthogonal and $\bar{Q}_{k} \bar{R}_{k}$ is the QR factorization of $A^{k+1}$ with $\bar{R}_{k}=R_{k} R_{k-1} \cdots R_{0}$.

  (b) Suppose that $A$ is symmetric and that its eigenvalues satisfy

  $$\left|\lambda_{1}\right|<\left|\lambda_{2}\right|<\cdots<\left|\lambda_{n-1}\right|=\left|\lambda_{n}\right|$$

  Suppose, in addition, that the first two canonical basis vectors are given by $\mathbf{e}_{1}=\sum_{i=1}^{n} b_{i} \mathbf{w}_{i}$, $\mathbf{e}_{2}=\sum_{i=1}^{n} c_{i} \mathbf{w}_{i}$, where $b_{i} \neq 0, c_{i} \neq 0$ for $i=1, \ldots, n$ and $\left\{\mathbf{w}_{i}\right\}_{i=1}^{n}$ are the normalised eigenvectors of $A$.

  Let $B_{k} \in \mathbb{R}^{2 \times 2}$ be the $2 \times 2$ upper left corner of $A_{k}$. Show that $d_{\mathrm{H}}\left(\operatorname{Spec}\left(B_{k}\right), S\right) \rightarrow 0$ as $k \rightarrow \infty$, where $S=\left\{\lambda_{n}\right\} \cup\left\{\lambda_{n-1}\right\}$ and $d_{\mathrm{H}}$ denotes the Hausdorff metric

  $$d_{\mathrm{H}}(X, Y)=\max \left\{\sup _{x \in X} \inf _{y \in Y}|x-y|, \sup _{y \in Y} \inf _{x \in X}|x-y|\right\}, \quad X, Y \subset \mathbb{C}$$

  [Hint: You may use the fact that for real symmetric matrices $U, V$ we have $\left.d_{\mathrm{H}}(\operatorname{Spec}(U), \operatorname{Spec}(V)) \leqslant\|U-V\|_{2} \cdot\right]$

  comment

• Paper 2, Section II, 41E

  Numerical Analysis

  Part II, 2021

  (a) Let $\mathbf{x} \in \mathbb{R}^{N}$ and define $\mathbf{y} \in \mathbb{R}^{2 N}$ by

  $$y_{n}= \begin{cases}x_{n}, & 0 \leqslant n \leqslant N-1 \\ x_{2 N-n-1}, & N \leqslant n \leqslant 2 N-1\end{cases}$$

  Let $\mathbf{Y} \in \mathbb{C}^{2 N}$ be defined as the discrete Fourier transform (DFT) of $\mathbf{y}$, i.e.

  $$Y_{k}=\sum_{n=0}^{2 N-1} y_{n} \omega_{2 N}^{n k}, \quad \omega_{2 N}=\exp (-\pi i / N), \quad 0 \leqslant k \leqslant 2 N-1$$

  Show that

  $$Y_{k}=2 \omega_{2 N}^{-k / 2} \sum_{n=0}^{N-1} x_{n} \cos \left[\frac{\pi}{N}\left(n+\frac{1}{2}\right) k\right], \quad 0 \leqslant k \leqslant 2 N-1$$

  (b) Define the discrete cosine transform $(\mathrm{DCT}) \mathcal{C}_{N}: \mathbb{R}^{N} \rightarrow \mathbb{R}^{N}$ by

  $$\mathbf{z}=\mathcal{C}_{N} \mathbf{x}, \text { where } z_{k}=\sum_{n=0}^{N-1} x_{n} \cos \left[\frac{\pi}{N}\left(n+\frac{1}{2}\right) k\right], \quad k=0, \ldots, N-1$$

  For $N=2^{p}$ with $p \in \mathbb{N}$, show that, similar to the Fast Fourier Transform (FFT), there exists an algorithm that computes the DCT of a vector of length $N$, where the number of multiplications required is bounded by $C N \log N$, where $C$ is some constant independent of $N$.

  [You may not assume that the FFT algorithm requires $\mathcal{O}(N \log N)$ multiplications to compute the DFT of a vector of length $N$. If you use this, you must prove it.]

  comment

• Paper 3, Section II, 40E

  Numerical Analysis

  Part II, 2021

  Consider discretisation of the diffusion equation

  $$\frac{\partial u}{\partial t}=\frac{\partial^{2} u}{\partial x^{2}}, \quad 0 \leqslant t \leqslant 1$$

  by the Crank-Nicholson method:

  $$u_{m}^{n+1}-\frac{1}{2} \mu\left(u_{m-1}^{n+1}-2 u_{m}^{n+1}+u_{m+1}^{n+1}\right)=u_{m}^{n}+\frac{1}{2} \mu\left(u_{m-1}^{n}-2 u_{m}^{n}+u_{m+1}^{n}\right), \quad n=0, \ldots, N,$$

  where $\mu=\frac{k}{h^{2}}$ is the Courant number, $h$ is the step size in the space discretisation, $k=\frac{1}{N+1}$ is the step size in the time discretisation, and $u_{m}^{n} \approx u(m h, n k)$, where $u(x, t)$ is the solution of $(*)$. The initial condition $u(x, 0)=u_{0}(x)$ is given.

  (a) Consider the Cauchy problem for $(*)$ on the whole line, $x \in \mathbb{R}$ (thus $m \in \mathbb{Z}$ ), and derive the formula for the amplification factor of the Crank-Nicholson method ( $\dagger$ ). Use the amplification factor to show that the Crank-Nicholson method is stable for the Cauchy problem for all $\mu>0$.

  [You may quote basic properties of the Fourier transform mentioned in lectures, but not the theorem on sufficient and necessary conditions on the amplification factor to have stability.]

  (b) Consider $(*)$ on the interval $0 \leqslant x \leqslant 1$ (thus $m=1, \ldots, M$ and $h=\frac{1}{M+1}$ ) with Dirichlet boundary conditions $u(0, t)=\phi_{0}(t)$ and $u(1, t)=\phi_{1}(t)$, for some sufficiently smooth functions $\phi_{0}$ and $\phi_{1}$. Show directly (without using the Lax equivalence theorem) that, given sufficient smoothness of $u$, the Crank-Nicholson method is convergent, for any $\mu>0$, in the norm defined by $\|\boldsymbol{\eta}\|_{2, h}=\left(h \sum_{m=1}^{M}\left|\eta_{m}\right|^{2}\right)^{1 / 2}$ for $\boldsymbol{\eta} \in \mathbb{R}^{M}$.

  [You may assume that the Trapezoidal method has local order 3 , and that the standard three-point centred discretisation of the second derivative (as used in the CrankNicholson method) has local order 2.]

  comment

• Paper 4 , Section II, 40E

  Numerical Analysis

  Part II, 2021

  (a) Show that if $A$ and $B$ are real matrices such that both $A$ and $A-B-B^{T}$ are symmetric positive definite, then the spectral radius of $H=-(A-B)^{-1} B$ is strictly less than $1 .$

  (b) Consider the Poisson equation $\nabla^{2} u=f$ (with zero Dirichlet boundary condition) on the unit square, where $f$ is some smooth function. Given $m \in \mathbb{N}$ and an equidistant grid on the unit square with stepsize $h=1 /(m+1)$, the standard five-point method is given by

  $$u_{i-1, j}+u_{i+1, j}+u_{i, j-1}+u_{i, j+1}-4 u_{i, j}=h^{2} f_{i, j}, \quad i, j=1, \ldots, m$$

  where $f_{i, j}=f(i h, j h)$ and $u_{0, j}=u_{m+1, j}=u_{i, 0}=u_{i, m+1}=0$. Equation $(*)$ can be written as a linear system $A x=b$, where $A \in \mathbb{R}^{m^{2} \times m^{2}}$ and $b \in \mathbb{R}^{m^{2}}$ both depend on the chosen ordering of the grid points.

  Use the result in part (a) to show that the Gauss-Seidel method converges for the linear system $A x=b$ described above, regardless of the choice of ordering of the grid points.

  [You may quote convergence results - based on the spectral radius of the iteration matrix - mentioned in the lecture notes.]

  comment